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I.  INTRODUCTION 

An  overall  rise  in  fuel  prices  has  led  to  an  increasing 
interest  in  the  design  of  autopilots  for  ships.  The  purpose 
of  the  automatic  steering  control  is  to  minimize  the  propul- 
sion losses,  which  are  caused  by  added  drag  due  to  steering 
of  the  ship.  minimizing  a  performance  criterion  based  on 
added  drag  due  to  steering  can  reduce  fuel  consumption. 
Claims  by  many  researchers  indicate  that  a  carefully 
designed  controller  could  save  from  one  to  two  percent  of 
fuel.  For  large  containerships  this  could  amount  to  more 
than  1100,000.00  per  year  savings. 

To  study  the  optimization  problem,  models  of  both  the 
ship  and  its  operating  environment  are  required.  What  type 
of  computer  model  should  be  used  to  represent  the  ship? 
Chapter  two  addresses  the  development  of  several  models. 
Since  the  best  model  was  desired  it  was  decided  to  use  the 
equations  of  motion  to  simulate  the  ship  in  our  Fortran 
program.  The  basic  Nomoto  models  give  an  adequate  descrip- 
tion of  ship  steering  dynamics  for  design.  The  Ncmoto 
second-  and  third-order  models  were  developed  from  the  equa- 
tions of  motion  as  defined  by  a  series  expansion  including 
all  terms  (both  linear  and  nonlinear)  for  which  hydrodynamic 
coefficients  were  available.  An  interactive  program  that 
utilized  the  Nomoto  mcdels  to  model  the  ship  was  also  used. 
Two  independent  programs  were  developed  to  aid  in  the  design 
of  the  controller. 

fihat  is  an  adequate  cost  function  which  represents  the 
added  drag  due  to  steering?  Chapter  three  addresses  the 
classical  cost  function  used  by  many  researchers. 

Since  a  variety  cf  control  algorithms  are  possible  one 
must  ask  if  one  algorithm  provides  a  lower  minimal  cost  than 
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another.  Chapters  fear,  five  and  six  address  the  selection 
of  the  controller  which  provides  the  minimum  value  of  added 
drag  due  to  steering. 

Ship  dynamics  change  with  operating  conditions  such  as 
ship  speed  ,  sea  state  ,  and  encounter  angle.  Therefore  an 
adaptive  controller  must  be  used  to  provide  minimum  added 
drag  due  to  steering.  Chapter  seven  development  of  an 
approach  to  an  adaptive  controller  utilizing  satellite 
information. 

Conclusions  were  drawn  from  these  experiments  and  are 
presented  in  Chapter  eight.  This  thesis  investigated  only 
course  keeping  with  emphasis  on  minimizing  rudder  and  yawing 
activity  to  reduce  fuel  consumption.  Presented  in  this 
Chapter  are  recommendations  for  future  study  where  the 
objective  is  track  fcllowing  which  would  be  important  for 
ships  reguired  to  follow  stringent  routes.  It  is  also  impor- 
tant for  other  systems  such  as  satellites,  missiles, 
aircraft,  where  the  controller  minimizes  yaw  error  to  keep 
the  system  on  track. 
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II.  COHPOTEB  MODELS 

The  aodel  which  test  represents  ship-steering  dynamics 
is  a  Taylor's  series  expansion  of  the  force  and  moment  rela- 
tionships  around  a  selected  steady-state  operating  point. 
The  resulting  equations  are  commonly  known  as  the  eguations 
of  motion  [Bef.  1].  A  computer  program  was  developed  using 
known  available  data  on  the  hydrodynamic  coefficients  for 
the  Sl-7  containership  to  provide  a  computer  simulation  of 
the  ship.  The  computer  program  is  shown  in  Appendix  A. 
figure  2.1  shows  the  block  diagram.  Small  yaw  command 
angles  are  used,  for  example  YAWC=  1.0  /  57.296  represents  a 
yaw  command  change  of  one  degree. 
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Figure  2.1    BLOCK  DIAGRAM 
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To  obtain  the  Noaoto  second-  and  third-order  transfer 
functions  from  the  equations  of  motion,  the  function  mini- 
mization subroutine  vas  used  to  obtain  the  coefficients. 
Figure  2.2  shows  the  scheme  used  to  obtain  the  Ncmoto 
transfer  functions-  The  computer  program  is  shown  in 
Appendix  A. 
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Figure  2.2    DETERMINATION  OF  NOMOTO  MODELS 

The  Soaoto  models  were  checked  against  analytic  results 
from  linearized  equations. 

Proceeding  to  the  second-order  Nomoto  eguation: 
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xP(S)/B(S)    =    K/S*(1+1*S)  (2.1) 

Deriving  the  second- ccder  Nomoto  transfer  function  from  the 
yaw    egaation   only,    the   result   is 

i(S)/8{S)    =    0.040893/S* (1+8.539932*S) 
and    using    function   mirimization   as   in   Figure   2-2 

+  [S)/b(S)     =    0.0409221/S*  (H8.5520782*S) 
and   the      agreement   is  obvious.      Osing      function   minimization 
with    both   yaw   and   sway    equations    with    linear   terms   only,    the 
results   are: 

^<S)/3(S)     =    0.  1072741/S*  (1+3  1.9  199524* S) 
If    the   nonlinear    terms   are      included   but   the   perturbation   is 
small 

+  {S)/6(S)     =    0.  1072082/S*  { 1 +3  1  .8907013*S) 
and    it    is   clear    that    the   nonlinear   terms   contribute    little. 

Proceeding    to   the   third-order   Nomoto  equation: 

i£(S)/5(S)    =    K*{1  +  TZ*S)/S*  (1+TP1*S)  *  (1+TP2*S)  (2.2) 

The  parameters  were  calculated  and  checked  by  using  function 
minimization  as  in  figure  2.2.  The  results  are  given  in 
Table  1.  It  is  clear  that  the  answers  obtained  by  function 
minimization   agree  closely   with   the   analytic   solutions. 


TABLE    1 
THIRD-ORDEB    NOHOTO    HODEL    FOR    THE    SL-7 

speed  K  T2  TP1  TP2 

knots   calc     comp     calc     comp     calc        comp          calc  comp 

16      .0738    .0738   22.57    22.95    12.946    12.946    107.583  107.583 

23       .1067    .1061    15.67    15.70      9.014      9.006      75.130  74.846 

32      .1477    .1477    11.28    11.28      6.470      6.467      53.793  53.793 
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Analytical      equations      used      to      calculate      seccnd-crder 
Nomotc   transfer   function    coefficients   are: 


K    =8      /N 

6         r 


T   =-  (I    -N.    )/N 
z       r  r 


Analytical      equations      used        to      calculate 
transfer   function   coefficients   are: 


third-crder 


K    =  (N    -N    *Y    /Y    )/(N    -H    *  {Y    -M*0)  /Y    ) 

6       v      6      vrvr  v 

TZ    =  -(  (E-Y-)  *N~  -N-*Y.  )/(Y    *N„-N    *ZX  ) 
*  *        v'       o       v      6  '     *    v      6       v      5  ' 

TP1*1P2=-  (  (M-Y-)*(I^-N-)-N-*Y-)/(N    *{Y    -M*U)-Y    *N    ) 
%*        v'     *  z      r'      v      r  v        r  v      r 

TPUTP2=(  {M-Y  •)  *N    ♦  (I   -N-)*Y    +  N.*{Y    -M*U)+Y-*N    ) 
*  v        vy      r     *  z      r        v      v    *   r  r      v 

/(N    *{Y    -M*U)-Y    *N    ) 
v  v    *   r  '      v      r 

The  nondimensionalized  hydrodynamic  coefficients  for  the 

SL-7  ccntainership  are  shown  in  Table  2. 


TABLE  2 
HYDBODYNABIC  COEFFICIENTS  FOB  THE  SL-7 


axial  force 


X  • 

X 

X 

x» 

X' 


uu 


vr 


vv 


66 


=  -0.0001 
-  -0,0003 
=  0.0039 
=  -0.0012 
=  -0.0005 


lateral  force 

Y«v   =  -0.00758 

V         =   0.0023 
r 

Y»    =   0.00145 

o 

Y«    =   0.01 
vvr 

Y»    =  -0.008 
vrr 

Y»    =  -0.03 
vvv 

Y'rrr=   0.003 
Y'665=  -0.0005 


moment  z-axis 

N»    =  -0.00213 

v 

N*    =  -0.00105 

r 

N»    =  -0.0007 

N'    =  -0.015 
vvr 

S»    =  -0.008 
vrr 

N»    =   0.01 
vvv 

N'    =  -0.006 
rrr 

H»    =   0.0001 

ooo 
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III.  COST  FUNCTION 

In  recent  years,  many  have  studied  the  problem  of 
[fief.  2]  [fief.  3]  [fief.  4]  [fief.  5]  [fief.  6]  [fief.  7] 
[fief.  8]  [fief.  9]  [fief.  10]  [fief.  11]  [fief-  12]  [Ref.  13] 
[fief.  14]  optimizing  an  automatic  ship-steering  controller 
for  minimum  fuel  consumption.  It  is  well  known  that  addi- 
tional drag  is  introduced  by  steering  and  that  both  the 
rudder  action  and  the  yawing  motion  contribute  to  this  added 
drag.  A  measure  of  the  added  drag  given  as  a  cost  function 
is 

T 
J    =    \/y  {  X*^2    +    52   )    dt  (3.1) 

o 

where   ^  =  yaw  error 

5  =  rudder  angle 

X  =  weighting  factor 

While  this  expression  is  an  approximation,  it  is  conven- 
ient for  shipboard  use  because  ^  and  5  are  readily  measur- 
able. There  is  no  general  agreement  on  numerical  values  for 
the  weighting  factor,  X  ,  and  in  this  study  the  values  used 
were  chosed  from  the  work  of  fi.E.  Reid  [fief.  7]  for  the 
SL-7. 

The  weighting  factors  for  the  operating  range  of  the 
ship  are  shown  in  Table  3. 
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TABLE  3 
WEIGHTING  FACTOH 

ship   speed  weighting   factor 

(knots) 

16  16.796 

23  8-128 

32  4.2 


Reid*s  work  shows  the  relationship  of  weighting  factor 
to  the  closed-loop  natural  frequency,  mass,  pivct  point, 
ship  speed,  x"Vr  ani*  x<  66  hy^rodynamic  coefficients.  It  is 
shown  in  Equation  3.2.  Reid  chose  a  closed-loop  natural 
frequency  of  0.05  rad/sec  which  experimentally  shewed  at 
this  frequency,  the  weighting  factor  in  the  cost  function, 
provided  good  representation  of  the  added  drag  due  to 
steering. 


X    =      2*ft*  (1  +  X»      )*(CP/L)*o>2    /(p/2)*(HXtt     *02)  (3.2) 

vr  dd 
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IV.    COHTBOLLEB    DESIGN    USING    ICSOS 

The  Interactive  Control  System  Optimization  and 
Simulation  (ICSOS)  package  finds  optimum  values  for  unknown 
(free)  parameters  in  a  control  system  design  problem  and/or 
performs  simulation  of  the  system.  An  example  of  usage  of 
ICSOS   is   shown   in   Appendix   B. 

In  preliminary  studies  ICSOS  was  used  with  Nomotc  models 
to  study  controller  characteristics  in  calm  water.  The  func- 
tion minimization  subroutine  adjusted  the  controller  parame- 
ters to  ninimize  the  cost  function.  Figure  4.1  shows  the 
scheme   used  to   evaluate  the   controller   parameters. 
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Figure    4.1         OPTIMIZATION    OF   CONTBOLLEH 
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Beid  £Bef.  7]  uses  the  second-order  Nomoto  model  of 
equation  2.1  for  the  SL-7  and  also  uses  a  controller 
described    by 

Gc  (£)     =    K1*(1+T1*S)    /     (1  +  T2*S)  (4.1) 

i 

His   results   are   given  in  Table   4. 

TABLE    4 
BEID'S    BESOLTS 

speed  plant  weighting  controller  gains 

knots  K  T  factor  K1  T1  T2 

16  0.1084    90.36       16.796         0,4556       89.51  10.06 

23  0.1556    64.67         8.128         0.3769      62.60         8.308 

32  0.2167    45.45        4.2  0.3188      44.92         7.066 

Using    this      plant   and   weighting      factor    values      but   applying 
ICSOS,    results    were    obtained   and   shown   on   Table   5. 

TABLE    5 
ICSOS    BESOLTS 

speed  plant        weighting  controller    gains  cost 

knots        K  T        factor  K1  T1  T2        J    min 

16         .1084    90.36     16.796    .454616    90.3459    10.0215    340.864 
23  .1556    64.67       8.128    -373171    64.6658       8.4640    139.9916 

32         .2167    45. 45       4.2         .318645    45.4475       7.0662      60.828 

In  each  case  the  controller  zero  (1/T1)  cancels  the 
plant  pcle  (1/T) .  Additional  experiments  consisting  of 
inserting  arbitrary  numbers  in  the  Nomoto  equation  and 
repeating  the  computer  run  indicated  that  this  will  always 
be  true.  That  is,  to  minimize  the  cost  the  plant  pole  is 
cancelled  and  a  new  pole  location  determined  with  appropri- 
ately  adjusted    gain. 
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The  siiiFle  controller  of  Eguation  4.1  is  an  arbitrarily 
chosen  structure.  To  determine  the  effects  of  more  complex 
controllers  three  additional  structures  were  chosen  as  shown 
in  Figure  4.2.  Each  of  these  was  used  with  the  Ncmoto 
second-order  model  for  the  ship  at  each  of  the  indicated 
speeds.  The  results  are  shown  in  Tables  6,  7 ,  and  8. 
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Figure  4.2    YlfilOOS  STBOCTOBES  FOB  COHTBOllEBS 

These  results  are  very  interesting.  At  16  knots  the 
controller  gain  (K1) ,  controller  zero  (1/T1)  and  controller 
pole  (1/12)  are  essentially  the  same  for  all  structures.  For 
structure  B,  which  includes  an  additional  pole,  the  function 
minimization  subroutine  tries  to  drive  the  additional  pole 
to  infinity,  and  no  doubt  would  have  done  so  if  the  calcula- 
tions had  continued.  For  structure  C,  which  has  two  poles 
and  two  zeros,  a  zero  and  pole  cancel  indicating  that  they 
are  not  needed  or  wanted.    For  structure  D,   the  integrator 
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TABLE  6 
SIflULATION  RESULTS  -  STEADY  STATE   600  SECS 

CALH  WATER  FOB  VABIOOS  CONTBOLLEBS 
FOB  FIXED  SHIP  SPEED  (  16  KNOTS  j 
HCHOTO  SECOND-ORDER  HODEL  (K=- 1084- T=90. 36) 
X=  16.796  ,    OPTIMAL  PARAMETER  GAINS  OF 
VABIOUS  CCNTBOLLEBS  ,  COST  FUNCTION 

contr  controller  gains  cost. 

K1      T1       T2      T3      T4     Ti     J  min 

A  .454616  90.435S  10.0215  -      -  340.864 

B  .444101  90.2950  9.8566  .01      -  -    341.046 

C  .454511  90.3685  10.0224  23.085  23.084  -    340.864 

D  .454581  90.371S  10.0222  -       -  1E09  340.864 


TABLE  7 
SIMULATION  RESULTS  -  STEAD!  STATE   600  SECS 

CALH  WATEB  FOB  VABIOOS  CONTROLLEBS 
FOB  FIXED  SHIP  SPEED  {  23  KNOTS  ) 
NGMOTO  SECOND-ORDER  MODEL  (K=. 1556 ,T=64. 67) 
X=  8.128  ,  OPTIMAL  PARAMETER  GAINS  OF 
VARIOUS  CCNTBOLLEBS  ,  COST  FUNCTION 

contr         controller  gains 

K1      T1         T2        T3       T4 

A   .373171  64.6657S   8.463957 

B   .340024  79.65872   8.889204    .01 

C   .373139  64.66855   8.463497  25.9719  25.9738 


cost 

J 

mm 

139. 
140. 
139. 

,9916 
,93  38 
.3991 

TABLE  8 
SIMULATION  RESULTS  -  STEADY  STATE   600  SECS 


CALH  WATEB  FOB  VABIOUS  CONTBOLLEBS 
FOB  FIXED  SHIP  SPEED 
NOMOTO  SECOND-OBDEB  MODEL 

X=  4.2    .  OPTIMAL  PARai 
VARIOUS  CONTROLLERS  .  COST  FUNCTION 


J32  KNOTS  J 
K=-2167.T=45.45) 
METER  GAINS  OF 


contr  controller  gains  cost 

K1  T1       T2        T3         T4  J  min 

A   .318645  45.44747  7.06617     -         -  60.828 

B   .318  45.45     7.066       .05        -  60.933 

C   .318678  45.57511  7.06790  50.1829  50.04832  60.828 
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gain  is  driven  to  zero.  The  same  pattern  of  results  is 
obtained  at  23  knots  and  32  knots-  Note  that  in  all  cases 
the  ninimum  cost  is  essentially  the  same,  as  would  be 
expected   since    all    controllers   are   the   same. 

Using  the  computer  method  of  Figure  4.1  and  the  Ncmoto 
third-order  models  of  Table  1 ,  controllers  A,  B,  C  of  Figure 
4.2  were  optimized..  The  results  are  shown  in  Tables  9,  10, 
and    11. 

TABLE    9 
SIMULATION    RESULTS   -    STEADI    STATE      600    SECS 

CALM    IATEB    FOB    VABIOUS    CONTBOLLEBS 
FOB    FIXED    SHIP    SPEED    f     16    KNOTS    ) 
NOMCTO    THIRD-ORDER    MODEL 
(K=. 0738 12, TZ=22. 5673, TP 1=1 2. 9458, TP2=107.  5853) 
X=    16.796    ,    OPTIMAL    PARAMETER    GAINS    OF 
VABIOUS    CCNTBOLLERS    ,    COST    FUNCTION 

contr  controller   gains  cost 

K1  T1  T2  T3  T4  J    min 

A       0.6446104       90.CS94       15.27712         -  -         370.4023 

B      0.6441367      84.826         15.78691    .24598  -  374.3808 

C      0.6151139    107.5782      8.73520    12.9368    24.9676    369.9297 


TABLE    10 
SIMULATION    RESULTS  -    STEADY    STATE      600   SECS 

CALM    WATER   FOR    VABIOUS    CONTBOLLEBS 
FOB    FIXEE   SHIP    SPEED    (    23   KNOTS   ) 
HOHCTO    THIRD-ORDER    MODEL 
(K=. 1067.TZ=15. 675.TP1=9.014, TP2=75.13) 
X=    8.128    ,    OPTIMAL    PABAMETEB    GAINS    OF 
VABIOUS    CCNTROLLERS    ,    COST    FUNCTION 

contr  controller  gains  cost 

K1  T1  T2  T3  T4  J    min 

A      0.5224258    63.13609    12.72212  -  -         152.2920 

B      0.5216467    64.93709    12.63218    .0505174  -         152.5333 

C       C. 5001907    75.14852    6.527490    9.039928    18.260    152.2800 


Of   uajcr     interest   is      the   fact      that    the     difference  in 
"cost"    between    A,      B,      C   is    less    than   one    per  cent.      At   each 
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TABLE  11 
SIMULATION  RESULTS  -  STEADY  STATE   600  SECS 

CALM  WATER  FOR  VABIOOS  CONTROLLERS 
FOR  FIXED  SHIP  SPEED  {  32  KNOTS  ) 
NOMOTO  THIRD-ORDER  MODEL 
(K=.  1477  1.TZ  =  11.  2833  ,TP1=6. 4699, TP 2=53.  7931) 
A.  =  4.2   -  OPTIMAL  PARAMETER  GAINS  OF 
VARIOUS  CONTROLLERS  ,  COST  FUNCTION 

contr  controller  gains  cost 

K1        T1        T2       T3        T4       J  min 

A  0.427633  48.66C48  10.74485  -  -  68.09039 
3  0.298732  89.40696  15.01033  .0597786  -  69.32355 
C   0.417991  53.69S61  4.970016  6.294354  13.85724  68.04735 


speed  (16,23,32  knots)  controller  C  is  "BEST",  bat  the 
difference  is  slight.  Examining  the  parameter  values 
obtained  for  controller  C,  it  is  seen  that  at  all  three 
speeds  both  poles  of  the  ship  are  essentially  cancelled  by 
zeros  of  the  controller. 

These  results  seem  to  indicate  that  the  dynamics  cf  the 
plant  determines  the  optimum  structure  for  the  controller. 

Using  a  state-feedback  controller  and  Nomoto  third-order 
models  of  Table  1,  the  controller  was  optimized  for  various 
ship  speeds.  Figure  4.3  shows  the  scheme  used  to  evaluate 
the  state-feedback  controller. 

Using  the  scheme  of  Figure  2.2,  with  no  change  in  cost 
function  or  weighting,  the  optimal  gains  and  costs  were 
deter lined  as  shown  in  Table  12. 

Hhen  comparing  the  state- feedback  controller  with 
controller  C,  it  is  seen  that  at  each  speed  controller  C  has 
a  lower  cost.  Among  the  controllers  consired,  controller  C 
is  "BEST"  when  using  the  Nomoto  third-order  model,  although 
the  differences  in  cost  are  not  dramatic. 
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TABLE  12 
SIBOLATION  BESOLTS  -  STEADY  STATE   600  SECS 

CAlfl  WATIB  POB  VABIOOS  SHIP  SPEEDS, 

OPTIBAL  P&BABETEB  GAINS  FOB 

STATE-FEEDBACK  CONTBOJLLEB 

speed   Hcmoto  third-crder  plant  weighting  controller  cost 
knots   K    TZ     TP1    TP2   factor    K1    K2   J  min 

16  .0738  22.567  12.946  107.583  16.796  4.426  78.004  370.711 
23  .1067  15.675  9.014  75.13  8.128  3.103  45.649  152.596 
32   .1477  11.283   6.470   53.793   4.2   2.240  27.896  68.2513 


25 


T.  COHTBOIIEB  DESIGN  OSING  FOBTBAN  PROGRAM 

The  Fortran  program  referenced  in  Chapter  two  which 
provided  a  computer  simulation  of  the  SL-7  ship  was  modi- 
fied. A  function  minimization  subroutine  was  coulped  to  the 
simulation  and  used  the  subroutine  to  adjust  controller 
parameters  to  minimize  the  cost  function  and  to  evaluate  the 
minimum  cost.  Figure  5.1  shows  the  scheme  used  to  evaluate 
the  controller  parameters.  This  prograa  was  used  for  compar- 
ison to  ICSOS.  The  computer  program  is  shown  in  Appendix  A. 
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Figure  5.1    OPTIMIZATION  OF  CONTROLLER  OSING  FORTRAN  PROGRAM 
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Using  the  computer  method  of  Figure  5. 1  and  the  nonli- 
near equations  of  motion,  controllers  A,  B,  C  of  Figure  4.2 
were  optimized.  The  results  are  shown  in  Tables  13,  14,  and 
15. 

TABLE  13 
SIBOLATIOH  RESULTS  -  STEADY  STATE   600  SECS 

CALM  WATER  FOR  VARIOUS  CONTROLLERS 
FOR  FIXED  SHIP  SPEED  (  16  KNOTS  ) 
EQUATIONS  OF  MOTION 
A=  16.796  ,  OPTIMAL  PARAMETER  GAINS  OF 
VARIOUS  CONTROLLERS  ,  COST  FUNCTION 

contr  controller  gains  cost 

K1      T1        T2         T3         T4       J  min 

A  .646401  89.81704  15.381699  -  -  1.128189 
B  -620050  90.67294  15.542297  0.9201336  -  1.173323 
C    .617326  107.1494   8.597198  13.353928  25.21362  1.126307 


TABLE  14 
SIMULATION  RESULTS  -  STEADY  STATE   600  SECS 

CALM  iATEH  FOR  VARIOUS  CONTROLLERS 
FOR  FIX1I  SHIP  SPEED  (  23  KNOTS  ) 
EQUATIONS  OF  MOTION 
X=   8.128  .  OPTIMAL  PARAMETER  GAINS  OF 
VARIOUS  CONTROLLERS  #  COST  FUNCTION 

contr  controller  gains  cost 

K1       T1       T2         T3       T4        J    min 

A    .522106  66.33122  12.83327  -       0.4640879 

B    .4S5869  66.15152  13.01183  0.92783      -       0.4857854 
C    .503967  74.79771   6.65880  9.20533  18.4022064  0.4636095 


These  results  agree  with  those  obtained  by  ICSOS  and 
controller  C  provides  the  minimum  cost. 

If  the  assumption  that  the  steering  dynamics  of  the  ship 
is  adequately  modeled  as  a  second-order  system  is  valid, 
then  only  two  states  are  needed  for  feedback.  For  a  third- 
order  system  three  states  are  required.  The  controller 
structures  are  shown  en  Figure  5.2  and  5.3.  Using  the  scheme 
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TABLE  15 
SIMULATION  RESULTS  -  STEADY  STATE   600  SECS 

CALB  BATES  FOB  VABIOUS  CONTROLLERS 
FOB  FIXES  SHIP  SPEED  (  32  KNOTS  ) 
EQUATIONS  OF  MOTION 
X=   4.2    .  OPTIMAL  PABAHETEB  GAIHS  OF 
VABIOOS  CCNTBOLLEBS  ,  COST  FUNCTION 


contr 


K1 


controller  gains 
T1       T2       T3 


T4 


cost 
J   fflin 


A       .428404    48.65540    10.814426  -  -  0.2072417 

B       .298732    89.40696    15.010330      0.01  -  0.2118334 

C       .417333    53.09654      5.096548    6.474857    14.0205    0.2071124 


of  Figure  5.1  with  no  change  in  cost  function  or  weighting, 
the  optimal  gains  and  cost  were  determined  as  shown  in 
Tables  16  and  17.  Coaparing  costs,  there  is  little  differ- 
ence between  the  twc  state  system  and  the  three  state 
system.  The        comparison      between     state        feedback      with 

controller  C,  it  is  seen  that  at  each  speed  controller  C  is 
tetter,    but  not    much   tetter. 
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Figure  5.3    THBEE  STATE  SYSTEH 


TABLE  16 
SIBOLATIOH  HISOLTS  -  STEADY  STATE   600  SECS 


speed 
knots 

16 
23 
32 


CALH  IATEB  FOB  VARIOUS  SHIP  SPEEDS 
EQUATIONS  OF  MOTION 
OPTIHAL  PARA  METER  GAINS  FOB 
TiO  STATE  SISTE3 


K1 


gains 


K2 


4.4033689 
3.0889006 
2.2342062 


77.5041656 
45.2637787 
27.6808014 


weighting 
factor 

16.796 
8.128 
4.2 


cost 
J   Bin 

1.128771 

.4646050 
.2075207 


Note  that  for  the  Nomoto  model  studies  yaw  error  and 
rudder  angles  vere  measured  in  degrees;  when  the  equations 
of  action  were  simulated  yaw  error  and  rudder  were  in 
radians.  Thus  the  numerical  values  of  the  cost,  J,  are 
different.  , 

Transient  response  plots  for  controllers  A,  B,  C,  and 
three  state-feedback  at  ship  speed  32  knots  are  shown  in 
Figures   5.4  -  5.9. 
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TABLE  17 
SIflOLATIOH  BESULTS  -  STEADY  STATE   600  SECS 


speed 
knots 

16 
23 
32 


CALH  BATEB  FOB  VARIOUS  SHIP  SPEEDS 
EQUATIONS  OP  HOTION 
OPTIBAL  PARAMETER  GAINS  FOB 
TBBEE  STATE  SYSTEM 


K1 


4.8617249 
3.6630983 
2.5967150 


gains 
K2 

87.7073364 
56.2784882 
33.7511444 


K3 

99.9802704 
88.5913391 
41.3186035 


weighting 
factor 

16.796 
8.128 
4.2 


cost 
J    sis 

1.12  82  89 

.46435  48 
.2074225 
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Figure   5.4        YAB   vs.    TIHE      (controller  A,    B,    C  and   state- feedback) 
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figure   5.5        RODDEB   *s»    TIME      (controller   A,    B,    C   and   state-feedback) 
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Figure   5.6        TAW    AND   HODDEE    vs.    TIME       (controller   A) 
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Figure   5.8        TA1   AHD   BUDDEB   vs.    TIBE       (controller  C) 
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Figure  5.9   YAW  AHD  BUDDEB  vs.  TIME   (state-feedback  controller) 
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VI.  CJWTJOLLEB  DESIGN  IN  SEA  STATE 

A  sea  state  generator  was  coupled  to  the  Fortran 
program,  so  that  the  function  minimization  subroutine  could 
be  used  in  the  presence  of  the  sea  state.  The  sea  state 
generator  was  an  elaborate  program  obtained  from  DTNSEDC. 
This  program  generates  added  mass  and  inertia  values  as 
functions  of  encounter  frequency  and  also  calculates  the 
forces  and  moments.  The  forces  and  moments  are  generated 
and  stored  in  a  look  up  table  which  was  coupled  to  the  equa- 
tions of  motion.  figure  6. 1  shows  the  scheme  used  to  eval- 
uate the  controller  parameters.  The  computer  prograi  is 
shown  in  Appendix  A. 

The  optimal  gains  obtained  by  the  calm  water  study  of 
Chapter  five  were  use  as  the  initial  guess  in  evaluating  the 
optimal  controller  parameters  in  the  presence  of  a  seaway. 
For  comparison,  studies  of  the  value  of  the  cost  function 
using  calm  water  gains  in  sea  state  were  obtained;  then  the 
function  sinimiza ticn  subroutine  was  allowed  to  adjust 
controller  parameters  in  the  presence  of  several  sea  states 
and  encounter  angles.  The  entire  study  was  done  at  a  ship 
speed  of  32  knots.  Ihe  added  mass  and  inertia  change  with 
respect  to  encounter  frequency  as  shown  in  Figures  6.2  and 
6.3.  Figure  6.3  is  ncndimensionalized  by  dividing  the  added 
inertia  by  the  mass  of  the  displaced  water  and  the  square  of 
the  length  between  perpendiculars.  To  convert  back  to  dimen- 
sionalized  units  of  lb-ft-sec2,  multiply  the  graph  points  by 
2. 581E12.  Since  the  sea  state  is  represented  by  irregular 
waves,  the  waves  impinging  on  the  ship  hull  contain  the 
total  energy  density  spectrum  composed  of  many  frequencies 
and  the  ship  responds  to  an  average  value  of  added  lass  and 
inertia.    The   values  used  for   this  study  was   obtained  at 
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encounter  frequency  cf  0.75  rad/sec  from  our  sea  state 
generator.  This  frequency  gave  us  values  for  added  mass  and 
inertia  representative  of  an  average  value.  The  energy 
density  spectra  for  various  sea  states  are  shown  in  Figure 
6.4.  The  added  mass  for  sway  was  changed  from  2.6457E06 
lb-sec2/ft  for  calm  water  to  2.3043E06  lb-sec2/ft  for  a 
seaway.  The  added  moment  of  inertia  for  yaw  was  changed 
from  1.42E11  lb-ft-sec2  for  calm  water  to  1.5096E11 
lb-ft-sec2  for  a  seaway.  All  other  hydrodynamic  coeffi- 
cients were  kept  constant  at  calm  water  values.  The  results 
are  shown  in  Tables  18  -  25.  In  certain  sea  states  and 
encounter  angles  the  calm  water  optimal  gains  performed  well 
as  shewn  by  calm  water  cost  value  when  compared  to  sea  state 
cost  value.  In  most  cases  the  function  minimization  suhrou- 
tine  found  new  gains  with  lower  cost  function  values  in 
seaway  as  compared  to  using  calm  water  gains.  In  the  calm 
water  evaluation,  the  system  was  perturbed  with  a  one  degree 
course  change,  but  the  course  change  was  not  included  in  the 
seaway  tests.  The  difference  in  cost  values  is  attributed  to 
the  difference  in  operating  conditions. 

Osing  the  Proportional,  Integral  and  Derivative  (PID) 
controller  Equation  6.1  with  no  change  in  cost  function  , 
the  function  minimization  subroutine  was  used  to  adjust 
controller  parameters  to  minimize  the  cost  function  and 
evaluate  the  minimum  cost.  The  results  are  shown  in  Table 
26.  fihen  comparing  the  PID  with  controller  A,  it  is  seen 
that  at  each  encounter  angle,  controller  A  is  better.  These 
results  agree  that  in  a  seaway  controller  A  provides  the 
minimum  cost. 

Table  27  shows  comparison  of  the  minimum  cost  function 
for  controller  A,  controller  C,  and  PID.  The  study  was  done 
at  ship  speed  of  32  knots  and  at  sea  state  4.  Controller  A 
provides  the  minimum  cost. 
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Pigure  6.1    OPTIMIZATION  OP  CONTBOXLEB  IN  SEA  STATE 

The  optimal  gains  obtained  in  the  presence  of  sea  state 
was  dene  over  using  a  simulating  time  of  600  seconds.  The 
sea  state  program  is  designed  to  provide  gradual  increase  in 
the  forces  and  moments  during  an  initial  time  interval.. 
This  is  done  to  minimize  initial  condition  transients  in  the 
ship  dynamics.  There  will- unavoidably  be  some  transient 
effects,  however,  and  these  could  affect  the  value  of  the 
cost,  J,  determined  during  the  600  seconds  of  simulation.  To 
determine  whether  such  initial  transients  had-  any  signifi- 
cant contribution  to  the  value  of  J,  additional  simulation 
runs  were  made  with  the  controller  parameters  fixed  at  their 
optimal  values.  However,  evaluation  of  the  cost,  J,  was 
started  only  after  300  seconds   of  simulation   had  elapsed. 
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Figure   6,2         AEDED    HASS    vs.    EHCOONTEB    FHEQOEHCI 
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TABLE  18 
SIMULATION  BESOLTS  -  STEADY  STATE   600  SECS 

FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  SEAWAY 
SHIP  MODEL:  EQUATIONS  OF  MOTION 
YAWC=0.0 

SEA  STATE  1 

CONTBOILEB    A 

encounter  controller   gains  sea   state        cost    with 

angle  cost  calm    water 

degree  K1  11  T2  J    min  J 

0  .4284037  48-6554395  10.814426  .617452-34  .61745E-34 

30  1.1561117  29.3693695  1.4592390  .2870198  .5128402 

60  1.4033298  10.6530075  1.1086683  .1342071  .2154726 

90  .2969198  58.2413940  1.8758221  .1300669  .1565958 

120  .1761794  299.999512  30.7967834  .05741726  .0727727 

150  2.8430557  5.2826872  .8887b96  .0219070  .0939400 

180  1.6211386  14.0782928  2.0712433  .0051925  .0095694 


TABLE    19 
SIBOLATION    HESOLTS   -    STEADY   STATE      600    SECS 

FIXED    SHIP    SPEED     (32    KNOTS)  IN    A    SEAWAY 
SHIP    MODEL:    EQUATIONS    OF    MOTION 
YAWC=0.0 

SEA    STATE   2 

CONTBOILEH    A 

encounter          controller   gains  sea   state        cost    with 

angle  cost               calm    water 

degree             K1                11                  T2  J   min                    J 

0         .42840370    48.6554395    10.814426  .61745E-34    .61745E-34 

30         .27997030    249.935059    19.857742  .04774852       .0886225 

60         .95575100    24.3813629    2.3079853  .04104504       .0535879 

90         1.3577642    9.43564080     1.1068363  .02650556       .0483197 

120         1.12C8973    25.4498596    4.0224676  .04928402      .0717524 

150          2.S777727    16.2154541     .56274800  7.5751530       28.1294403 

180         .61420630    .482041200    6.2521963  .000124338    .0002445 
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TABLE  20 
SIMULATION  BESULTS  -  STEADY  STATE   600  SECS 

FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  SEAWAY 
SHIP  MODEL:  EQUATIONS  OF  MOTION 
YAWC=0.0 

SEA  STATE  4 

CONTBOLLEB    A 

encounter  controller   gains  sea  state        cost    with 

angle  cost  calm    water 

degree  K1  11  12  J    min  J 

0  .4284037  48.65540  10.814426    .620598E-34  .620598E-34 

30  .9815440  5.7^3036  .6999879       .02854677  .0395892 

60  .6201209  40.8C556  19.606873    .09375697  .1032696 

90  1.809746  36.C1225  6.324703       1.5171340  4.1623011 

120  5.195190  18.92513  .6999907       9.991730  48.970703 

150  1.446776  16.89375  .5265408       16.67052  24.822098 

180  .1000000  1.000000  20.149999    .00739631  .0076657 


TABLE    21 
SIMULATION    BESOLTS   -    STEADY    STATE      600   S2CS 

FIXED    SHIP    SPEED    (32    KNOTS)     IN    A    SEAWAY 
SHIP    MODEL:    EQUATIONS    OF    MOTION 
YAWC=0.0 

SEA    STATE   6 

CONTBOLLEB    A 

encounter  controller  gains  sea   state        cost   with 

angle  cost  calm   water 

degree        K1  11  T2  J   min  J 

0  .4284037      48.6553955    10.8144264  .50399E-32    .50899e-32 

30  2.9715786    10.4721832.5342450  1.4287940    4.74724010 

60  1.7228041     8.4014740       .5141125  1.5827220    3-42744920 

90  1.6584366    37.1672655    .5792384  4.5505371     13.2757149 

120  3.3422489    106.722259    .9260592  22.103002    94.5497589 

150  .2854474       157.483887    119.981018  .81100580     1.50448510 

180  .8053379       .75733550      6.04484460  .07365978    .142564400 
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TABLE    22 
SIMULATION    BFSULTS   -    STEAD!    STATE      600    SJ5CS 

FIXED    SHIP    SPEED    {32    KNOTS)     IN    A    SEAWAY 
SHIP    MODEL:    EQUATIONS    OF    MOTION 
YAWC=0.0 

SEA    STATE    1 

CONTBOILER  C 

encounter  controller   gains  sea  calm 

angle  cost  cost 

degree      K1  T1  T2  T3  T4  J    min  J 

0  1.61345    16.8755  25.0481  47.3405  1.73759    .  14E-33  .45E-33 

30  .957558    12.6178  7.32113  43.3531  8.15752    .324739  .357733 

60  .781984    17.6475  9.22485  13.9438  16.7663    .159710  .203137 

90  .417332    53.0965  5.09655  6.47857  14.0205    .148588  .148588 

120  .417332    53.0965  5.09655  6.47857  14.0205    .077918  .077918 

150  2-13735    18.8265  17.5778  25.1516  21.1481    .031496  .081612 

180  .957558    12.6178  7.32113  43.3531  8.15752    .006566  .008172 


TABLE    23 
SIMULATION    BFSULTS   -   STEADY    STATE      600    SECS 

FIXED    SHIP    SPEED    (32    KNOTS)    IN    A    SEAWAY 
SHIP    MODEL:    EQUATIONS    OF   MOTION 
YAWC=0.0 

SEA    STATE   2 

CONTBOILEB    C 

encounter  controller   gains  sea  calm 

angle  cost  cost 

degree      K1  T1  T2  T3  T4  J    min  J 

0  .417333    53.0965  5.09655  6.47436  14.0205    .18E-33  . 1 8E-33 

30  .849594    19.9913  9.34138  20.3578  13.7487    .054338  .061184 

60  .417333    53.0965  5.09655  6.47436  14.0205    .044536  .044536 

90  .781984    17.6475  9-22485  13.9438  16.9438    .033518  .048467 

120  .880395   21.4597  10.9255  9.24547  11.0667    .055292  .056288 

150  .899999    15.7103  1.11632  41.3275  2.29012    10.7636  23.8522 

180  .440916    .093671  17.7305  25-2103  5.04178    .000125  .001635 
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TABLE  24 
SIMULATION  BESULTS  -  STEADY  STATE   600  SECS 

FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  SEAWAY 
SHIP  HODEL:  EQUATIONS  OF  MOTION 
YAHC=0.0 

SEA  STATE  4 

CONTBOLLEB    C 

encounter  controller   gains  sea  calm 

angle  cost  cost 

degree      K1  T1  T2  T3  T4  J    min  J 

0  1.22424   70.3578  61.3016  10.5467  61.8215    . 62E-34    .142-33 

30  .690573    20.1488  20.3214  5.10369  19.7841    .034033    .071978 

60  .782547    12.6178  13.7713  21.5637  21.5637    .098914    .244369 

90  2.22895    51.7744  54.3190  17.3522  6.07814    1.57368    2.98305 

120  3.72749   85.4697  40.6999  8.52234  1.35207    10.3530    37.3988 

150  .417333    53.0965  5.09655  6.47486  14.0205    20.3956    20.3956 

180  .059166    .286208  19.5103  46.4767  22.3286    .007397    .099413 


TABLE   25 
SIMULATION    BESULTS   -    STEADY    STATE      600    SECS 

FIXED    SHIP    SPEED    (32    KNOTS)    IN    A    SEAHAY 
SHIP    MODEL:    EQUATIONS    OF    MOTION 
YAHC=0.0 

SEA    STATE   6 

CONTBOLLEB    C 

encounter  controller   gains  sea  calm 

angle  cost  cost 

degree      K1  T1  T2  T3  T4  J    min  J 

0  2.33178   52.0881  95.7892  31.6564  11.6959    .49E-32  .19E-31 

30  2.08709   73.1270  76.7193  12.1726  16.4711    2.00375  3.83379 

60  2.00128    71.6612  77.6750  13.1170  17.0912    2.15428  3.41794 

90  .957558    12.6178  13.7713  72*0670  15.4611    5.76399  8.80971 

120  3.10589    81.8044  38.2439  91.2237  9.14683    24.9099  72.6716 

150  1.51250   70.3578  61.3016  35.5894  61.8215    2.50971  7.50022 

180  1.52875    1.87828  43.4698  49.9147  11.2599    .078894  .930885 
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TABLE  26 
SIHDLATION  BESOLTS  -  STEADY  STATE  600  SECS 

FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  SEAWAY 
SHIP  MCEEL:  EQDATIONS  OF  MOTION 
YAiC=0-0 

SEA  STATE  4 

PID   CONTROLLER 

encounter  controller   gains  sea  state 

angle  K1  T1  T2  J    min 

30  .95263100  4.20720860  .69368610  .02965619 

60  .68631890  12.5794449  8.2121658  .09730512 

90  2.5809155  12.4247589  .77810380  1.5915950 

120  4.9198265  12.5986176  .67592390  10.708980 

150  1.3970823  15.7682953  .51991180  17.4272C0 


&(S)/+    (S)    =    K1    +    K1*T1*S    /     (UT2*S}**2  (6.1) 


TABLE  27 
COMPARISON  OF  TEE  MINIMUM  COST 


SHIP  SPEED  (32  KNOTS) 
SEA  STATE  4 


encconter 
angle 

controller 
A 

controller 
C 

controller 
PID 

degree 

J  min 

J  min 

J  min 

30 

60 

90 

120 

150 

.02854677 
.093756S7 
1-5171340 
9.99173CC 
16.67052G 

.034033 
.098914 
1. 57368 
10.3530 
20.3956 

.029  85619 
.09730512 
1.5915950 
10.708980 
17.427200 

The  value  obtained  was  then  doubled  and  compared  with  the 
result  of  evaluating  J  over  the  full  600  seconds.  Comparison 
of  Table  28  with  cost  values  in  Tables  18,  19,  20,  21  shows 
only  small  differences. 

To  obtain  insight  into  the   stochastic  process  of  irreg- 
ular seas,  a  deterministic  process  was  studied.   The  Fortran 
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TABLE  28 
EFFECTS  DDE  TO  TRANSIENT  AND  GRADUAL  BUILD  OP  OF  SEA  STATE 

INTEGRATION  OF  COST  FUNCTION  (  300  TO  600  SECS) 
FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  SEAWAY 

sea  encounter      controller  gains        cost     ccst 
state  angle     K1        T1        T2      J  min    2*J  min 

1  60  1. 4033296  10.650075  1.1086683  .0641122  .12-82244 

2  60  .955751C0  24.381363  2.3079853  .0199731  .0399462 
4  60  .62012090  40.805560  19.606873  .0515974  .1031948 
6  60  1.7228041  8.4014740  .51411250  -7906  179  1.581236 


program  was  modified  to  minimize  the  cost  function  in  the 
presence  of  a  regular  sea.  To  allow  comparison  with 
previous  work  the  encounter  frequency  of  0.75  rad/sec  was 
used  and  scaled  the  amplitude  of  the  regular  sea  to  its 
prospective  sea  state.  The  entire  study  was  done  at  a  ship 
speed  of  32  knots.  The  results  are  shown  in  Ta-bles  29/  30, 
and  2  1. 

Table  29  shows  that  for  regular  seas  the  ccntrcller 
parameters  do  not  change  significantly  for  different  sea 
states;  tut  as  sea  state  increases,  the  cost  value  increases 
due  to  the  increase  in  yaw  moment  and  sway  force  on  the 
ship.  Tables  30  and  31  also  show  that  the  controller  parame- 
ters do  not  change  significantly  from  sea  state  to  sea 
state.  However,  an  encounter  angle  of  90  degrees  shews  a 
relatively  high  cost  compared  to  costs  calculated  for  6C  and 
120  degrees  at  a  given  sea  state.  To  account  for  this 
anomaly,  the  following  is  suggested.  In  the  regular  sea, 
the  added  mass  and  inertia  were  known  for  a  given  encounter 
frequency,  while  in  the  irregular  sea  a  representive  average 
value  was  used.  The  nethod  used  to  obtain  the  average  might 
not  represent  the  actual  average.  Also,  it  seems  reasonable 
to  suppose,  that  the  assumptions  of  the  function  weighting 
factor  are  satisfied  for  all  encounter  angles;  that  is,  the 
weighting  function   (Eg.   3.2),   which   appears  in   the  cost 
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function  {Eq-  3.1),  does  totally  represent  the  added  drag 
for  all  encounter  angles.  Future  study  is  needed  tc  answer 
these  questions. 

The  sea  state  in  the  deterministic  model  is  represented 
by  regular  waves.  On  this  description,  the  waves  impinging 
on  the  ship  hull  correspond  to  only  one  frequency  in  the 
energy  density  spectrum.  In  the  case  of  irregular  seas, 
however,  the  spectral  components  change  for  different 
states,  as  shown  in  Figure  6.4.  Thus  comparison  of  the 
controller  parameters  obtained  for  regular  seas  with  results 
for  irregular  seas  is  not  justified.  The  function  miniiiza- 
tion  subroutine  adjusted  controller  parameters  to  minimize 
the  cost  function  fcr  either  case  (irregular  or  regular 
seas)  as  shown  in  tables  32  and  33. 

TABLE  29 
SIHULATION  RESULTS  -  STEADY  STATE   600  SECS 

FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  REGULAR  SEAWAY 
SHIP  HGDEL:  EQUATIONS  OF  MOTION 
YAHC=0.0 
ENCOONTEE  FREQUENCY  =  0.75  RAD/SEC 
ENCOUBTER  ANGLE  =  60  DEGREES 

CONTROLLER  A 

controller  gains  cost 

T1  T2  J  min 

141.383179  32.9405670  .000764582 

129.987473  31.4042358  .003056434 

135.798737  32.9749756  .009345479 

135.488495  33.5585632  .022174600 


sea 
state 

K1 

1 
2 

4 
6 

.1449795 
.1534657 
.1514665 
.1533340 

Note  that  in  both  the  deterministic  and  stochastic 
models,  among  the  controllers  considered,  controller  A  is 
"BEST"  in  a  seaway  disturbance,  although  the  differences  in 
cost  are  not  dramatic. 

Finally,  the  observed  dependence  of  optimal  controller 
gains   on  sea   state  and  encounter  angle   suggests  that  an 
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TABLE  30 
SIMULATION  RESULTS  -  STEADY  STATE   600  SECS 

FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  REGULAR  SEAWAY 
SHIP  flCDEL:  EQUATIONS  OF  MOTION 
YAHC=0.0 
ENCOUHTEB  FREQUENCY  =  0.75  RAD/SEC 
SEA  STATE  4 

CONTROLLER  A 

encounter  controller  gains  cost 

angle       K1  T1  T2  J  min 

30  .2351043  102.021973  28.3396912  .002985300 

60  .1514665  135-798737  32.9749756  .009345479 

90  .4964442  66.546493  49.7598267  .048143090 

120  .1327230  149.540543  33.6013489  .038937880 

150  .4536914  70.566528  31.5839539  .062534153 


TABLE  31 
SIMULATION  RESULTS  -  STEADY  STATE   600  SECS 

FIXED  SHIP  SPEED  (32  KNOTS)  IN  A  REGULAR  SEAWAY 
SHIP  MODEL:  EQUATIONS  OF  MOTION 
YAWC=0.0 
ENCOUHTEB  FREQUENCY  =  0.75  RAD/SEC 
SEA  STATE  6 

CONTROLLER  A 

encounter  controller  gains  cost 

angle       K1  T1  T2  J  min 

30  .2370022  100.122940  28.0581207  .007092878 

60  .1533340  135.488495  33.5535632  .022174600 

90  .5210407  62.153702  49.9358093  .112772880 

120  .1414837  142.695160  35.3171234  .091541650 

150  .4587426  71.451385  33.4568024  .144615829 


adaptive  controller   nust  be   used  to   provide  a   continuous 
minimum  on  the  cost  function. 

After  obtaining  tre  optimal  gains  for  controller  A,  to 
observe  the  behavior  of  the  rudder  and  yaw  motion  of  the 
ship,  transient  resccnse  plots  were  obtained  for  controller 
A  at  ship  speed  of  32  knots  and  sea  state  4  for  various 
encounter  angles  as   shown  in  Figures  6.5  -   6.14.   Note  the 


49 


TABLE  32 
CCHPABISCN  OF  IRREGULAR  TO  EEGDLAE  SEAS  CONTROLLER  GAINS 

SEA  STATE  4 

CONTROLLER    A 

encounter  controller    gains 

angle                                                  K1                          T1  T2 

30              (irregular)         .9815440             5.733036  .6999879 

30                 (regular)            .2351043         102.021973  28.3396912 

60              (irregular)         .6201209           40.805560  19.6068730 

60'              (regular)            .1514665         135.798737  32. 974S756 

90              (irregular)         1.809746           36.012250  6.3247080 

90                 (regular)            .4964442           66.546493  49.7598267 

120              (irregular)         5.195190           13.925130  .6999907 

120                 (regular)            .1327230         149.540543  33.6013489 

150              (irregular)         1.446776            16.893750  .52654C3C0 

150                 (regular)            .4536914           70.566528  31.5839539 


TABLE    33 
COHPARISCN   OF    IRREGULAR    TO    REGULAR    SEAS    CONTROLLER    GAINS 


SEA    STATE 

6 

CONTROLLER 

A 

encounter 
angle 

controller    gains 
K1                          T1 

T2 

30 
30 

(irregular) 
(regular) 

2.9715786 
.23700220 

10-4721832 
100. 122940 

.534  2450 
28.05812C7 

60 
60 

(irregular) 
(regular) 

1.7228041 
.  15333400 

8.4014740 
135.488495 

.5141125 
33.5565632 

90 
90 

(irregular) 
(regular) 

1.8584366 
.5210407 

37-1672655 
62.153702 

.5792384 
49.9858093 

120 
120 

(irregular) 
(regular) 

3.3422489 
.  1414837 

106.722259 
142.695160 

.9260592 
35.3171234 

150 
150 

(irregular) 
(regular) 

.2854474 
.4587426 

157.483887 
71.451385 

119.981018 
33.4568024 

increase  in  both  rudder  and  yaw  amplitude  as  the  encounter 
angle  increased.  This  is  due  to  the  increase  in  yaw  moment 
and   sway  force    on    the  ship. 
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32  KNOTS-SEfl  STATE  4 
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Figure  6.5    YAH  vs.  TIME    30  DEGREES 
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Figure  6.6    BODDEB  vs.  TIME    30  DEGBEES 
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32  KNUTS-SCfl  STATE  4 


ENCOUNTER  ANGLE  60  DEGREES 
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Figure  6.7       IAI  vs.   TIME 
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Figure  6.8   RUDDER  vs.  TIME   60  DEGREES 
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Pigure  6.9    TAB  vs.  TIME 
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Figure  6.10    RODDER  VS.  TIME    90  DEGREES 
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32  KNOTS-SEfl  STATE  4 
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Pigure   6.11        XAI    vs.    TIHE         120    DEGREES 
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32  KNOTS-SEA  STATE  4 
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Figure   6. 12        EODDEH    vs.    TIME      120    DEGREES 
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Figure  6.13    YAH  vs.  TXUE    150  DEGBEES 
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Figure  6.14    BODDEE  vs.  TIME   150  DEGREES 
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VII-  AN  ADAPTIVE  CONTROLLER 

In  a  seaway,  the  controller  gains  changed  dramatically 
for  changes  in  sea  state  and  encounter  angle.  An  adaptive 
controller  must  be  used  to  provide  continuous  operation  on  a 
minimum  cf  the  cost  function.  This  Chapter  addresses  a 
theoretical  design  cf  an  adaptive  controller. 

In  the  future,  there  will  be  tetter  measurement  cf  navi- 
gation than  can  be  provided  by  conventional  equipment  on 
board  a  ship*  Presently  the  Navy  is  involved  in  a  prcgraa 
that  will  provide  precision  navigation  data.  The 
NAVSTAR/GICBAL  POSITION  SYSTEM  (GPS)  £Ref.  15]  [Ref.  16] 
[Ref.  17]  will  provide  extremely  accurate  three-dimensional 
position  and  velocity  information  to  users  anywhere  in  the 
world.  The  position  determinations  are  based  on  the  measure- 
ment of  the  transit  time  of  RE  signals  from  four  satellites 
of  a  total  constellation  of  eighteen.  This  system  is  sched- 
uled to  be  fully  operational  in  1988.  At  present  (1984) 
there  are  four  NAVS1AR/GPS  satellites  in  operation  which 
allows  three  to  four  hours  per  day  of  navigation  tiaie. 
Already  the  Texas  Instrument  Company  markets  a  receiver  for 
this  system  where  GPS  can  be  used. 

The  Navy  Remote  Ocean  Sensing  System  (NROSS)  [Ref.  18] 
will  be  able  to  determine  wind  velocities  over  the  world's 
oceans  with  an  accuracy  sufficient  to  determine  ccean 
surface  waves.  It's  objective  will  be  to  acquire  global 
ocean  data  for  operation  and  research  use  by  both  the  lili- 
tary  and  civil  sectors.  This  system  is  scheduled  to  launch 
its  first  satellite  in  June  1989. 

The  scheme  for  an  adaptive  controller  is  shown  in  Figure 
7.1.  Having  stored  the  optimal  controller  parameters  in  a 
look  up  table   as  functions  of  ship  speed,    sea  state,   and 
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encounter  angle,  the  ship  operating  condition  must  be  known 
so  that  the  table  is  useful.  NAVSTAB/GPS  would  identify 
ship  speed  and  NEOSS  would  identify  sea  state  and  encounter 
angle.  The  optimal  parameters  can  then  be  looked  up  and 
inserted  into  the  controller.  This  should  place  system  oper- 
ation near  the  minimun  J.  To  ensure  fine  tuning,  a  micro- 
processor programmed,  with  the  function  minimization  or-line 
in  machine  language,  with  inputs  of  yaw  error  and  rudder 
motion  of  the  ship  would  accomplish  the  fine  tuning  rapidly. 
Since  the  subroutine  is  written  in  Fortran  (as  used  for  this 
study)  this  would  be  inappropriate  for  on-line  use. 

The  adaptive  controller  can  be  performed  with  digital 
circuits  rather  than  analog  components.  Garcia  [Eef.  19] 
demonstrates  the  process  for  converting  an  analog  controller 
into  a  digital  controller.  Figure  7.2  illustrates  the 
processing  of  the  major  components  in  a  digital  controller. 
An  analog  component  circuit  can  be  replaced  by  an  analog  to 
digital  converter,  a  digital  processor,  and  a  digital  to 
analog  converter.  Soire  of  the  benefits  which  can  be  realized 
by  doing  this  are: 

1.  A  high-speed  processor  could  actually  process  a 
number  of  multiplexed  signals,  performing  processing  func- 
tions on  a  number  of  independent  channels. 

2.  The  processing  function  is  permanent  in  software, 
unless  deliberately  changed,  and  will  not  drift  with  age. 

3.  The  processing  function  can  be  changed  without 
changing  components,  merely  by  changing  software. 

4.  Accuracy  can  be  made  very  high  and  can  be  changed 
merely  by  changing  software. 

5-  Processing,  which  previously  required  large  compo- 
nents such  as  inductors  in  low-f reguency  controllers,  can 
now  be  performed  by  very  small  digital  circuits. 
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Figure   7, 1         ADAPTIVE    CONTROLLER 
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Figure  7.2    DIGITAL  BLOCK  DIAGBAH 

In  converting  an  analog  controller  to  a  digital 
controller,  the  process  can  be  broken  down  into  the 
following  steps: 

1.  Determine  the  desired  analog  transfer  function. 

2.  Set  the  sampling  frequency. 

3-  Apply  the  bilinear  z-transformation. 

4.  Match  one  point  in  the  s  domain  to  the  z  domain. 

5.  Obtain  the  optimum  constant  coefficients. 

6.  Obtain  the  digital  transfer  function. 

7.  Ottain  the  simulation  diagram. 

The  optimal  controller  parameters  can  be  stored  in 
memory.  Intel  company  markets  a  4  megabit  non-volatile 
read/ write     bubble      memory.        It     is     supported      by     a     YSLI 
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controller  which  provides  a  black  box  interface.  It  is  easy 
to  use  and  can  be  used  with  any  8-  or  16-bit  microproces- 
sors,  Tte  bubble  memory  advantage  is: 

1.  Fast  access  time  compared  with  disk  or  tape, 

2-  Ncn-volatile. 

3.  aide  temperature  range  of  operation, 

» 

4.  forking  storage. 

5.  Portable  operation 

6.  Lew  power, 

7.  High  reliability. 
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VIII-  CONCLUSIONS  AND  RECOMMENDATIONS  FOB  FOTOBE  STOCI 

A-   CONCLUSIONS 

In  designing  the  controller,  three  different  ship  ncdels 
were  used.  Using  the  second-order  Nomoto  model  Equation  2.1 
allowed  comparison  of  results  with  Reid's  [Bef.  7]  £Bef-  10] 
work.  It  is  clear  that  the  answers  obtained  by  function 
linimization  agree  closely  with  field's  results  as  shewn  in 
Tables  4  and  5-  A  tetter  description  of  the  ship  is  the 
third-order  Nomoto  model  which  involves  both  the  sway  and 
yaw  eguations.  This  nodel  includes  the  two  dominating  poles 
of  the  ship.  The  best  model  to  describe  the  dynamics  of  the 
ship  is  a  Taylor's  series  expansion.  This  allowes  both 
linear  and  nonlinear  terms  in  the  eguations  of  motion  to 
affect  tie  design  of  the  controller. 

To  determine  which  controller  structure  would  provide 
the  minimum  cost  due  to  steering,  various  structures  were 
studied.  It  was  found  that  the  dynamics  of  the  plant  deter- 
mines the  optimum  structure  for  the  controller.  In  calm 
water  study,  when  using  a  second-order  Nomoto  model,  the 
best  structure  was  controller  A.  When  the  third-order  Nomoto 
model  Eguation  2.2  was  used  the  best  structure  was 
controller  C,  but  the  difference  is  slight.  Observe  that  in 
each  case  the  controller  zeros  cancel  the  plant  poles.  When 
the  eguations  of  motion  were  used  for  the  plant,  the  best 
structure  was  controller  C.  When  the  eguations  of  motion 
were  coupled  to  a  sea  state  generator  and  the  cost  function 
was  ttinittized  in  the  presence  of  a  seaway,  the  best  struc- 
ture was  controller  A.  This  study  concludes  that  controller 
A  should  be  used. 

A  function  minimization  subroutine  is  an  engineer's  tool 
which  can  be  used  in  nany  engineering  problems.  Previously  a 
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ssatchcd  filter  was  designed  for  the  Naval  Postgraduate 
School  research  project  on  the  Space  Transportation  System 
(STS)  for  the  Get  Away  Special  Program.  It  was  matched  to 
the  signature  of  the  auxiliary  power  unit  (APU)  on  hoard  the 
space  shuttle.  The  goal  was  to  turn  on  the  solid  state 
recording  system  before  lift  off,  to  record  the  acoustic 
power  generated  inside  the  shuttle  bay.  Basically  the 
matched  filter  is  a  finite  Impulse  Response  (FIE)  filter 
with  the  weights  calculated  to  obtain  the  least  squared 
error  of  the  desired  output  when  the  input  is  the  signature 
of  the  APO.  Figure  8.1  shows  the  scheme  used  to  evaluate 
the  FIR  weights. 

B.   EICCHMENDAIIONS  FCR  FOTOEE  STUDY 

In  the  future  most  ships  both  military  and  commercial 
will  have  GPS  receivers  as  part  of  their  navigation  equip- 
ment. Using  extremely  accurate  three-dimensional  position 
and  velocity  information  from  satellite  platforms  will  allow 
ships  to  navigate  accurately  in  and  out  of  ports.  The  func- 
tion ainimization  subroutine  is  a  powerful  tool  for 
designing  the  controller.  This  routine  simply  takes  the 
inputs  that  require  minimization  and  adjusts  the  parameters 
to  accomplish  this  task.  The  cost  function  for  the  added 
drag  due  to  steering  is  a  function  of  yaw  error  and  rudder 
motion.  The  use  of  function  minimization  and  NAVSTAE/GPS 
provides  the  means  for  optimization  for  guidance  and  cont- 
roll.  There  are  several  areas  that  need  future  study  and 
work. 

1.  Should  the  objective  change  to  track  following  then 
it  is  necessary  to  minimize  the  yaw  error  only.  This  would 
be  very  important  both  militarily  and  commercially  should  a 
port  be  mined.  If  the  ship  could  follow  a  stringent  route, 
knowledge  of  mine  locations  would  allow  access. 
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Figure  8.1    HATCHED  FILTER  DESIGN 

2.  A  controller  for  orbit  keeping  for  satellites  with 
High-Energy  Laser  weapons  would  be  very  important.  The  small 
far- field  spot  size  cf  a  focused  laser  beam  can  be  selec- 
tively focused  on  the  most  vulnerable  component  on  the 
target,  Facilitating  precision  energy  deposition,  and 
greatly  increasing  the  probability  of  a  kill. 

3.  An  adaptive  controller  to  minimize  track  error  on 
board  a  cruise  missile  could  be  programmed  for  selective 
targets. 
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4.   Military  and  commercial  aircraft  can  benefit  just  as 
do  ships  ty  reducing  drag  to  minimize  fuel  consumption- 


/ 


69 


APPENDIX    A 
PROGRAM    10    CALCULATE    OPTIMAL    GAIHS 

The  program  is  set  up  to  calculate  the  optimal  gains  for 
controller  A.  It  is  referenced  in  Chapter  five  and  six.  It 
can  easily  be  modified  to  obtain  optimal  gains  for  the  rest 
of  the  controllers.  After  obtaining  the  optimal  gains  the 
program  most  be  modified  to  do  a  simulation.  The  program  has 
sufficient  comments  for  appropriate  changes.  It  is  refer- 
enced   in   Chapter    two. 

This      program      car.      be   modified      to      obtain      the      Ncmcto 
models.      It   is    referenced    in   Chapter    two.    The  following   need 
to   be   changed. 
C       GAIN    COEFFICIENTS    TO    BE    OPTIMIZED 
K=XX{1) 
TP1=XX  (2) 
TZ=XX(3) 
TP2=XX  (4) 
C       EBRCB    SIGNAL    TO    EEIVE    RUDDER     (YAW    ACTUAL    -    YAW    CCMMANE) 
C      FOR    EQUATIONS    OF    KGTION. 

D=YAW    -    YAWC 
C       ERROR    SIGNAL    TO    IBIVE   RUDDER     {YAW    COMMAND   -    YAH    ACTUAL) 
C       FCB    KCMOTO    3  ED    OEIEfi    MODEL. 
E2=YAWC-YAW2 
X1=  (D2-X2)/TF1 
X3=K*(TZ*X1  +  X2) 
X4=(X3-X5)/TE2 
C       INTEGRATION 

X2=X2+X1*DELI 
X5=X5+X4*DELT 
YAW2=YAW2+X5*EELT 
C      CCST    FUNCTION 

TDIFF=TDIFF+     (YAW-YAW2) **2 

A 
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PROGRAM  TO  CALCULATE  OPTIMAL  GAINS  FOB  CONTBOILEB 


JOB  (2  220,0356) , » RESEARCH ■ ,CLASS=J 

ON=1024K 


//GARCIA 

//  EXEC  FRTXCLGP,IMSI=DP,REGI 

//FCBI.SYSIN  DD  * 
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TBIS  PROGRAM  Will  OBTAIN  T 
IT  IS  EEFERENCED  IN  CHAPTE 
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CHANGE  XS(*)  TO  X( 
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XS 
XI 

xu 


Pi 

:i 


XU  ( 
1)=.427   '  '  ' 
2)=48.66 
3} =10.7 
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THE  LONER  LIMIT  F 
THE  UPPER  LIMIT  F 
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xs 

xs 
xs 
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IS 
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XL 
XU 
XL 
XU 
XL 
XU 
A  DESCEIFT 
IS  DISCUSSED 
E=9./13. 
NTA=1000 
NPE=100 
NAV=0 
NV=3 
IP=0 
TEE  FOLLOWING 
CAIL  SLANT (X) 

IE    ONLY    SIMULATION    IS    WANT 
CALL    BQXPLX(NV,NAV,NPR 
WBITS     (6,25) 
FORMATM*,1     OPTIMAL    GA 
DC    30    1=1.3 
HBITE(6. 40)  I,3S(I) 
FORM  AT  (ix/x  (,#I2,*)=,  ,F14.  7) 
STOP 
END 

SUBROUTINE  P 
SUBROUTINE  PLANT 
COMMON  TDIFF 


OR  THE  I'TH  VARIABLE 
OR  THE  I'TH  VARIABLE 


STATEMENT  MUST  BE  CHANGED  TO 


ED 
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LANT  (XX) 
(XX)  SIMULATES  THE  SHIP 


REAL*8 
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REAL*8 
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EEAL*8 
REAL*8 
DIMENSION' XX'(3) 


L,L2,L3,I4,L5,L6 
X,XDOT,S, YDOT.U, 
TIME,ETIME,XUDOT 
YV,YR,YE,YVVR, YV 
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INITIALIZE    THE    COST    FUNCTION 

ISE=0.0 

ISR=0.0 

TDIFF=0.0 

LAMDA=4.2 


EGRATION 
ONDS 
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c 
c 


c 
c 
c 


c 
c 


c 
c 
c 
c 


D    = 


ANGLE 


GAIN    COEFFICIENT 

K1=XX  (1 

T1=XX  [2 

T2=XX(3 
X,XDCI,Y,Y 

X  =  0-0 

Y  =  0.0 
XECT=0. 0 
YDCI=0.0 

U,UDCT,V,VDOT    Afi 

V  =  0.0 
UEOT=0.0 
VDCT=0. 0 
XA«=O.0 
E=0.0 
REOI=0.0 

CBDEEEE    SPEED    IN 
54. 01    FT/SEC=32 
UC=54.01 
AT    STEADY    STATE 

u=uc 

RUDDER 

D=0.0 

1=880.5 

L2=L**2 

L3=L*L*L 

L4=L*L3 

L5=L*L4 

L6=L*L5 
SEA  DISTURBANCE 
FCRCES  IN  X,Y  DI 
MOMENTS  IN  2 

FX=0. 

FY=0. 

MZ=0. 
ISEA    IS    A    SWITCH 

ISEA=1 
HYEBCDYNAMIC    COEF 

EBC=1.9876 

MASS=(.  0044)  * 

12=  (0.00028)  * 

YAWE=0.0 

X2=0.0 

DX2=0.0 

CONTINUE 

S=DSQET (U**2+ 
INPUT  YAW  COMMAN 

YAWC=0.0 

IF  (TIME.GE.O 
EBBCE  SIGNAL  TO 
(  COMPENSATOR  FI 

YAWE=YAW    -    YA 

DX2=(YAWE-X2) 

D  =  K1*(T1*DX2  + 
AXIAL  FORCE  HYDR 
XUDOT  IS  THE  ADD 
DIFFERENT    ENCOUN 

XUDOT=(-.000  1 
XU=  (-0.0253)  * 
XDU=(-0.0003) 
XVR=  (0.0039)  * 
2VV=    -.0012)* 
XDD=(-0.0005) 
LATERAL    FORCE    HY 
YV=  (-0.00758) 
YB=  (0.0023)  *( 
YE=  (0.00145) * 
YVVE=(0.01)  *( 


S    TO    BE    OPTIMIZED 


OT    ARE    FIXED    COORDINATES    ON    EARTH 


E    FIXED    COORDINATES    ON    SHIP 


FEET/SEC 
KNOTS 

ACTUAL    SPEED     (U)     =    COMMAND    SPEED     (UC) 


RECTICN    COMPUTED    IN    POUND    FORCE 

;    ISEA=0     (CALM    WATER)     ISEA=1      (SEA   STATE) 

FICIENTS    ARE    INSERTED    HERE    AS    PARAMETERS 

.5*RHO*L3) 
.5*RHO*L5) 


200 


V**2) 

D 


.0) 
DEI 

LIE 
WC 

X2) 

CEY 
EE 

1ER 


YAWC=0.0 
VE    RUDDER (YAH    ACTUAL    -    YAM    ORDERED) 

B    ) 


NAMIC    COEFFICIENTS     (SURGE) 
MASS    TERM    WHICH    MUST    BE    CHANGED    FOE 
ANGLES    ,    SPEED    ,    ENCOUNTER    FREQUENCY 


)*(.5*RHO*L3) 
(.5*RHO*L2*S) 
*(.5*RHO*L2) 
(.5*RHO*L3) 
(.5*RHO*L2$ 

*  (.5*RHO*L2*S**2) 
EfiODYNAMIC    COEFFICIENTS     (SSAY) 

*  J.  5*RHG*L2*S) 
,5*RHO*L3*S) 
(.5*RHO*L2*S**2) 
.5*RHO*L3/S) 
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YVRR=  (-0.008)  *  (.  5*RHO *L4/S) 

YVVV=    -0.03)  *  (.5*RHO*L2/S) 

IHJBE=  0.003)*    .5*RHO*L5/S) 

YDDD=  (-0-0005)*{-5*EHO*L2*S**2) 
C       YUDOT    IS    THE    ADDEE    MASS    TEBM    WHICH    MOST    BE   CHANGED    FOR 
C       DIFFERENT    ENCOUNTER    ANGLES,    SPEED    ,     ENCOUNTER    FREQUENCY 
C 

C  SVDOT=(-0.0039)*  (.5*RHO*L3) 

C    SPEED=32    KNOTS. ENCOUNTER    ANGLE    =150    , ENCOUNTER    FREQ    =.75 

XVDOT=-2.3043E+C6 
C       MCMENT    ABOUT    Z-AXIS    HYDRODYNAMIC    COEFFICIENTS     (YAW) 

NV=(-0.  00213)  *{.5*RHO*L3*S) 

NB=  (-0.00  105    *  t.5*RHO*L4*S) 

NB=  (-0.0007)  *  (.5*RHO*L3*s**2) 

SVVR=(-0.015)  *{.5*RHO*L4/S) 

NVRE=    -0.008) *]. 5*RHO*L5/S) 

NVVV=    0.0  1l*(.5*RHO*L3/S) 

NRBE=    -0.006)  *  l.5*RHO*L6/S) 

NDDD=jO.  0  001)_*    . 5*RHO*L3*S**2) 
C    NRDOT    IS    THE    ADDEE   INERTIA    TERM    WHICH    MUST    BE   CHANGED    FOR 
C    DIFFERENT    ENCOUNTER    ANGLE    ,    SPEED    ,    ENCOUNTER    FREQUENCY 
C 

C  NRDOT=(-0. 00027)  *(.5*RH0*L5) 

C       SEEEC=32    KNOTS. ENCOUNTER    ANGLE    =150    , ENCOUNTER    FREQ    =.75 

NBDOT=-1.5096E+11 
C       SETS    SEA    STATE    TO   ZERO 

IF     (ISEA.EQ. 1)    GO    TO    30 

FX=0. 

FY=0. 

MZ=0. 

GC  TO  35 
C   TABLE  LOOK  UP  OF  £EA  DISTURBANCE, 
C   UKIT  12  HAS  THE  SEA  STATE  LATA  NAMED  CH 
C   IT  MUST  BE  SYNCHRCNIZED  BY  APPROPRIATELY 
C   CALLING  CH  IN  THE  PROPER  TIME  IN  THE  LOOP. 
C   TEE  SEA  DATA  WAS  CREATED  FOR  600  SECONDS 
C   WITH  AN  INCREMENTAL  INTERVAL  OF  1  SECOND. 
30    READ  (12)  CH 

FX=CH  (3" 

FY=CH(4 

MZ=CH  (8' 
35         CCNTINU! 
C       U    ACTUAL   SPEED 
C       UC    COMMANDED    SPEEE 
C       XE    =    EEOPELLER    THEUST 

XE=-XUU*UC**2 
C      EQUATIONS    OF    MOTION 
C       FOB    CONSTANT    SPEEE    COMMENT    THE    NEXT    TWO   INSTRUCTIONS 

UDOT=(     (XVR    +    £ASS)*V*R    +    XUU*U**2    ♦    XVV*V**2 
1    ♦    XDD*D*D    +       EX    ♦    XP    ) /(MASS-XUDOT) 

VDOT=(YV*V    +     nR-MASS*U)*R    ♦    YD*D    +    YVVR*V**2*R 

1  +    YVRR*V*R**2    +    YVVV*V**3 

2  ♦    YRRR*R**3    +    YDDD*D**3    ♦    FY    ) /  (MA SS-YVDOT) 
BEOT=  (    NV*V    +    KR*B    «■    ND*D    ♦    NVVR*V**2*R 

1  +    NVBR*y*R**2    +    NVVV*V**3 

2  +    NRRR*B**3    +    NDDD*D**3    ♦    MZ    )/(IZ-NRDOT) 
C       WEEN    TO    PRINTOUT 

IF     (ICOUNT.EQ.11)     GO    TO    50 
GC   TO    300 
C      CCNVEET    RADIANS    TC    DEGREES 
50       YAWDEG=    YAW*57.296 
RDEG=R*57.296 
RDDEG=RDOT*57.296 
DDEG=D*57.296 
YAKC=YAWC*57.296 
C  WRITE    (6,100)     1IME,XP,X,XD0T.Y,YD0T 

C  1       ,UC,U€UDOT, V,VDOT,YAWC,YAWDEG,RDEG,RDDEG,EDEG 

100    FCRMAT-llX,,TIME=,,F8.3f  »    SEC       XP=',FlO-2,»     LBF       X= 
1#F8.2,*    FT      XECT=» ,F8.4,»     FT/SEC      Y=',F8-2, 
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2»    FT      YDOT=»,F8.4,»    FT/SEC  ,/,2x.  '     UC='    F8.4, 
31    FT/SEC         U=»,F8.4,'     FT/SEC         UDOT= ' , FlO .6 , 
4    f    FI/SEC**2         V^.FS.a,'     FT/SEC  VDOT=  ■  ,  F  10-  6  , 

5»    FT/SEC**2» ,/,2X* YAWC=» ,F8.4,»     DEG       YAW  =  ',F15.7, 
6'    DEG         YAW    RATE=*    F15.7,  ■     DEG/SEC         YAW    ACCEL=» 
7,F15.7,«     DEG/£EC**2',/,2X, 'RUDDER    =  «,F15.7,«    DEG    •  ,/) 
ICCUNT=1 
C       TEST    IF    WANT    TO    STOP 

300       IF     (TIME.GE.  ETIME)     GO    TO    400 
C       INTEGRATION    STEP    SIZE    DE1T 

DELT=1.0 
C       INTEGRATION 

G=U+UDOT*DELT 
V=V+VDOT*DELT 
R=B+RDOT*DELT 
YAW=YAW+R*DEIT 
X2=X2+DX2*DELT 
C      CCNVERT    SHIP    TO    11X1^5   COORDINATES    ON    EARTH 
XDOT=U*DCOS  (YAK)  -V*DSIN  (SAW) 
YDOT=U*DSIN  (YAW)  +V*DCOS  (YAW) 
X=X+XDOT*DELT 
Y=Y+YDOT*DELT 
IIME=TIME+DELT 
ICCUNT=ICOUNT+1 
ISE=ISE    +    LAMDA*YAWE**2 
ISR=ISR    +    D**2 
GG    TO    200 
C  J=    TDIFF    =    COST    FUNCTION 

400       IEIFF=ISE-HSR 

WRITE J6, 5 00)     ISE,ISR,TDIFF,K1,T1,T2 
500      FCRHAT(*     » . 1 X  ,  'ISZ='     F1 5- 7, '       ISR='     F15.7    »       TOTAL=l 
1,Fl5.7,2x, *K1=»,F15.7,2X,fT1=»,F15.7, 2X,,T2=* ,F15.  7) 
REWIND    12 
RETURN 
EKD 
C       DELETE    ALL    THE    FCILOWING    SUBROUTINE    IF    SIMULATION    ONLY 
C       AND    NOT    OPTIMIZATION    IS    WANTED 
C . 

c 

C  SCEROUTINE    BOXELX  (CATEGORY    HO) 

C 

C  PURPOSE 

c 

C  BCXPLX  IS  A  SUBROCTINE  USED  TO  SOLVE  THE  PROBLEM  CF 

C  lccacting  A  MINIMUM  (OR  MAXIMUM)  OF  AN  ARBITRARY  OBJECT- 

C  iV€  function  SUBJECT  TO  ARBITEARY  EXPLICIT  AND/CR 

C  implicit  constraints  by  tHE  COMPLEX  METHOD  OF  M.J.  BCX. 

C  explicit  constraints  are  dEFINED  AS  UPPER  AND  LOWER 

C  bounds  on  the  independent  variables  IMPLICIT  constraints 

C  may  be  arbitrary  function  of  the  varlABLES.   TWO  FUN- 

C  ction  subprogram  tc  evaluate  the  objective  FUNCTION  AND 

C  implicit  conSTRAINTS.  RESPECTIVELY,  must  be  SUPPLIED 

C  by  the  user  (see  EXAMPLE  BELOW).  BOXPLX  ALSO  HAS  tHE 

C  option  to  perform  integer  programming,  where  the  values 

c  of  the  independent  variables  are  restricted  to  integers. 

C 

C      USAGE 

C 

C         CALL  BOXPLX  (NV, NAV, NPR ,NTA,R, XS, IP, XU ,XL, YMN  ,IER) 

C 

C  DESCRIPTION  OF  PARAMETERS 
C 

C  NV    AN  INTEGER  INEUT  DEFINING  THE  NUMBER  OF  INDEPENDENT 

C  VARIABLES  OF  THE  OBJECTIVE  FUNCTION  TO  BE  MINIMIZED. 
C    NOTE:   MAXIMUM  NV  +  NAV  IS  PRESENTLY  50.   MAXIMIM  NV  IS 

C  25.   IF  THESE  LIMITS  MUST  BE  EXCEEDED,  PUNCH  A  SOURCE 

C  DECK  IN  THE  USUAL  MANNER,  AND  CHANGE  THE  DIMENSION 

C  STATEMENTS. 
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c 

C  NAV  AN  INTEGER  INEUT  DEFINING  THE  NOHBES  OF  AUXILIARY  var 

C  iaELES  THE  USER  WISHES  TO  DEFINE  FOR  HIS  OWN  CONVENIENCE. 

C  TYPICALLY  HE  MAY  WISH  TO  DEFINE  THE  VALUE  OF  EACH  IMPLICI 

C  CONSTRAINT  FUNCTION  AS  AN  AUXILIARY  VARIABLE.   IF  THIS 

C  IS  DCNE,  THE  OPTIONAL  OUTPUT  FEATURE  OF  BOXPLX  CAN  BE 

C  USED  TO  OBSERVE  THE  VALUES  OF  THOSE  CONSTRAINTS  AS  THE 

C  SCIUTICN  PROGRESSES.   AUXILIARY  VARIABLES,  IF  USED, 

C  SHCULD  BE  EVALUATFE  IN  FUNCTION  KE  (DEFINED  BELOW). 

C  NAV  MAY  BE  ZERO. 

C 

C  NPR   INPUT  INTEGER  CONTROLLING  THE  FREQUENCY  OF  OUTPUT 

c  desired  for  diagnostic  purposes. 

C  IF  NPR  .IE.  0,  NO  CUTPUT  WILL  3E 

C  PROEUCED  BY  BOXPLX.   OTHERWISE,  THE  CURRENT  COMPLEX  CF 

C  K=  2*NV  VERTICES  AND  THEIR  CENTROID  WILL  BE  OUTPUT  AFTER 

C  EACH  NPR  PSRMISSIBIE  TRIALS.   THE  NUM3ER  OF  TOTAL  TRIALS, 

C  NUMEEE  OF  FEASIBLE  TRIALS,  NUMBER  OF  FUNCTION  EVALUATIONS 

C  AND  NUMBER  OF  IMPLICIT  CONSTRAINT  EVALUATIONS  ARE  IN- 

C  CIUDED  IN  THE  OUTPUT. 

C  ADDITIONALLY,  (WHEN  NPR  .GT.  0)  THE  SAME  INFORMATION 

C  WILL  BE  CUTPUT: 

C 

C  1)  IF  THE  INITIAL  EOINT  IS  NOT  FEASIELE, 

C  2   AFTER  THE  FIRST  COMPLETE  COMPLEX  IS  GENERATED, 

C  3   IF  A  FEASIELE  VERTEX  CANNOT  BE  FOUND  AT  SOME  TRIAL, 

C  4}  IF  THE  OBJECTIVE  VALUE  OF  A  VERTEX  CANNOT  BE  MADE 

C     NC-LCNGEE-WORST. 

C  5)  IF  THE  LIMIT  ON  TRIALS  (NTA)  IS  REACHED  AND, 

C  6)  WEEN  THE  OBJECTIVE  FUNCTION  HAS  BEEN  UNCHANGED  FOR 

C  2*NV  TRIALS,  INDICATING  A  LOCAL  MINIMUM  HAS  BEEN 

C  FOUND. 

C 

C  IF  TEE  USER  WISHES  TO  TRACE  THE  PROGRESS  OF  A  SOLUTION, 

C  A  CHOICE  OF  NPR  =  25,  50  OR  100  IS  RECOMMENDED. 

C 

C  NTA   INTEGER  INPUT  OF  LIMIT  ON  THE  NUMBER  OF  TRIALS 

c  allowed  in  the  calculation. 

C  IF  THE  USER  INPUTS  NTA  .LE.    0,  A  default 

C  VAIUE  CF  2000  IS  USED.   WHEN  THIS  LIMIT  IS  REACHED 

C  CGNTFCL  RETURNS  TO  THE  CALLING  PROGRAM  WITH  THE  BEST 

C  ATTAINED  OBJECTIVE  FUNCTION  VALUE  IN  YMN,  AND  THE  BEST 

C    ATTAINED    SOLUTION    POINT    IN    XS. 

C 

C  R    A  REAL  NUMBER  INPUT  TO  DEFINE  THE  FIRST  RANDOM  NDMEER 

C  USED  IN  DEVELOPING  THE  INITIAL  COMPLEX  OF  2*NV  VERTICIES. 

C  (0.  .GT.  R  .LT.  1.1  IF  R  IS  NOT  WITHIN  THESE  BOUNDS, 

C  IT  WILL  BE  REPLACE!  BY  1./3.  . 

C 

C  XS    INPUT  REAL  ARRAY  DIMENSIONED  AT  LEAST  NV+NAV. 

c  the  first  nv  must  contain  a 

C  FEASIBLE  ORIGIN  FCB  STARTING  THE  CAL- 

C  CULATICN.   THE  LAST  NAV  NEED  NOT  BE  INITIALIZED.   UECN 

C  RETURN  FROM  BOXPLX.  THE  FIRST  NV  ELEMENTS  OF  THE  ARRAY 

C  CONTAIN  THE  COORDINATES  OF  THE  MINIMUM  OBJECTIVE 

C  function,  AND  THE  REMAINING  NAV  (NAV  .GE.  0)     CONTAIN  THe 

C  values  OX  THE  CORRESPONDING  AUXILIARY  VARIABLES. 

C 

C  IP    INTEGEE  INPUT  FOR  OPTIONAL  INTEGER  PROGRAMMING. 

C  if  ip=1,  THE  VALUES  OF  THE  INDEPENDENT  VARIABLES  WILL 

C  be    replaced  WITH  INTEGER  VALUES  (STILL  STGRED  AS  REAI*4) . 

C 

C  XU    A  REAL  ARRAY  DIMENSIONED  AT  LEAST  NV  INPUTTING  THE 

C  upper  BOUND  ON  EACH  INDEPENDENT  VARIABLE,  (EACH  EXPLICIT 

C  CONSTRAINT).  INPUT  VALUES  ARE  SLIGHTLY  ALTERED  BY  BOXPLX. 

C 

C    XL         A    REAL    ARRAY    DIMENSIONED    AT    LEAST    NV    INPUTTING    THE 

c   lower   bound    on   each   independent 

C    VARIABLE,     (EACH    EXPLICIT    CONstraint) . 
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C  NOTE:   FOR  BOTH  XU  AND  XL  CHOOSE  REASONABLE 

C  VALUES  I?  NONE  ARE  GIVEN.  NOT  VALUES  WHICH  ARE 

C  magnitudes  ABOVE  CR  BELOW  THE  EXPECTED  SOLUTION. 

C  input  values  are  SLIGHTLY  ALTERED  BY  BOXPLX. 

C 

C    YMN      THIS    OUTPUT    IS    THE    VALUE     <REAL*4)     OF    THE    OBJECTIVE 

C    funcTICN,COERESPONEING    TO    THE    SOLUTION    POINT    OUTPUT    IN    XS 

C 

C  IER   INTEGER  ERROR  RETURN-   TO  BE  INTERROGATED  UPCN 

C  return  FROM  BOXPLX.   IER  WILL  BE  ONE  OF  THE  FOLLOWING: 

C 

C  =-1   CANNOT  FIND  FEASIBLE  VERTEX  OR  FEASIBLE  CENTROID 

C  AT  TBE  START  OR  A  RESTART  {SEE  'METHOD1  BELOW). 

C  =0    FUNCTION  VALUE  UNCHANGED  FOR  'N*  TRIALS.    (WHERE 

C  N=6*NV+10)   THIS  IS  THE  NORMAL  RETURN  PARAMETER. 

C  =1    CANNOT  DEVELCE  FEASIBLE  VERTEX. 

C  =2    CANNOT  DEVELOP  A  NO-LCNGER-WORST  VERTEX. 

C  =3    LIMIT  ON  TRIALS  REACHED.    (NTA  EXCEEDED) 

C  NOTE:   VALID  RESULTS  MAY  BE  RETURNED  IN  ANY  OF  THE 

C         ABOVE  CASES. 

C 

C   EXAMPLE  OF  USAGE 

C 

C  THIS  EXAMPLE  MINIMIZES  THE  OBJECTIVE  FUNCTION  SHOWN  IN 

C  the  EXTERNAL  FUNCTION  FE(X).   THERE  ARE  TWO  INDEPENDENT 

C  varlAELES  X(1)  &    X  (2)  ,  AND  TWO  IMPLICIT  CONSTRAINT 

C  function  X(3)  S  XJ4)  WHICH  ARE  EVALUATED  AS  AUXILIARY 

C  variarles  (see  EXTERNAL  FUNCTION  KE(X)  ). 

C 

C    DIKEKSICN      XS  (4)  ,XU  (2)  ,XL  (2) 

CC 

CC       S1AETING    GUESS 

C  XS 

C  X 

CC       UPP 

C  XU  (1)    =    6.0 

c  XU(2i    =    6.0 

CC       lOWEc    LIMITS 

C  XIM)    =    0.0 

C  XL   2)    =    0.0 

CC 

C  R  =    9./13. 

C  NTA    =   5000 

C  NPR    =    50 

C  NAV    =   2 

C  NV   =    2 

C  IP    =    0 

CC 

C  CALL    BOXPLX     (NV, NAV, NPR .NTA, R, IS, IP, XU, XL, YMN, IER) 

C  WRITE(6,1)      HXS(I)  ,1=1,4)  ,  YMN,  IER) 

C  1FCRMAT     [////,  •         THE    POINT    IS    LOCATED    AT     (XS  (I) =)     • 

c     2.  4(e13.7,5x)  ,// , 

C     3»    AND  THE  FUNCTION  VALUE  IS  »,E13.7,»  IER  =  ',15) 

CC 

C  STOP 

C  END 

CC 

CC 

C  FUNCTION       KEfX) 

C    EVALUATE    CONSTRAINTS.       SET    KE=0    IF    NO    IMPLICIT    CONSTRAINT 

C    is    viCLATED«    OR    SET    KE=1     IF    ANY    IMPLICIT 

c   constraint    is   violated. 

C  DIMENSION    X  (4) 

C  X1    =    X(1) 

C  X2   =    X{2) 

C  KE   =    0 

C  X(3)     =    X1    ♦     1.732051*X2 

C  IF     (X(3j     .LT.     0.     .OR.     X  (3)     .  GT.    6.)     GO   TO    1 

C  X  (4)    =    X1/1. 732051    -X2 
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:s(i)   =   1.0 
;s  2)    =   0.5 

'ER    LIMITS 
XUM)    = 
XU(2)    = 

iwec  ti a: 


C  IF  (X(4)  .GE.  C.)   RETURN 

CC 

C  1  KI  =  1 

C  RETURN 

C  END 

CC 

cc 

C      FUNCTION   FE(X) 
C      DIMENSION  X(4) 
CC 

CC   THIS  IS  THE  OBJECTIVE  FUNCTION. 

C      FE=  -{X(2)**3  *{9.-(X  (1)-3.)  **2)/(46.  76538)) 

C      RETURN 

C      END 

C 

C      METHOD 

C 

C  THE  COMPLEX  METHOE  IS  AN  EXTENSION  AND  ADAPTION  OF 

c  tJbe  simple  method  cf  linear  programming. 

C  STARTING  WITH  ANY  CNE  feasible  point  in  n-dimensicn 

C  A  "COMPLEX"  OF  2*N  vertices  is  constructed  by 

C  SEIECTING  RANDOM  ECINTS  WITHIN  THE  feasible 

C  REGICN.   FOR  THIS  PURPOSE  N  COORDINATES  ARE  FIRST 

C  RANDOMLY  CHOSEN  WITHIN  THE  SPACE  BOUNDED  BY  EXPLICIT  CCH- 

C  STRAINTS.   THIS  DEFINES  A  TRIAL  INITIAL  VERTEX. 

c  it  is  then  checked  for  possible  violation 

C  OE  IMPLICIT  CONSTRAINTS.   IF  one  or  more  are  violated, 

C  THE  TRIAL  INITIAL  VERTEX  IS  DISPLACED  half  of  its 

C  DISTANCE  FROM  THE  CENTROID  OF  PREVIOUSLY  SELECTED  initial 

C  VERTICES.   IF  NECESSARY  THIS  DISPLACEMENT  PROCESS  IS 

C  REPEATED  UNTIL  THE  VERTEX  HAS  BECOME  FEASIBLE.   IE  THIS 

c  FAIL   TO  happen  after  5*n+10  displacements, 

C  THE  SOLUTION  IS  AEANDONED.  AFTER  EACH  VERTEX  IS  ADDED 

C  TC  THE  COMPLEX,  THE  CURRENT  centroid  is  checked  fcr 

C  FEASIBILITY.   IF  II  IS  INFEASIBLE,  the  last  trail 

C  VERTEX  IS  ABANDONEE  AND  AN  EFFORT  TO  GENERATE  an  alter- 

C  ATIVE  TRIAL  VERTEX  IS  MADE.   IF  5*N+10  VERTICES  ARE 

C  ABANDCNED  CONSECUTIVELY,  THE  SOLUTION  IS  TERMINATED. 

C 

C  IF  AN  INITIAL  COMPIEX  IS  ESTABLISHED,  THE  BASIC 

c  computation  loop  is  initiated. 

C  THESE  INSTRUCTIONS  FIND  THE  CURRENT  WORST  vertex,  that 

C  IS,  THE  VERTEX  WITH  THE  LARGEST  CORRESPONDING  value  for 

C  THE  OEJECTIVE  FUNCIIGN,  AND  REPLACE  THAT  VERTEX  BY 

C  ITS  OVIR-REFLECTICN  THROUGH  THE  CENTROID  OF  ALL  CTHER 

c  vertices,  (if  the  vertex  to  be 

C  REPLACED  IS  CONSIDERED  AS  A  VECTOR  IN  n-space, 

C  ITS  OVER-REFLECTION  IS  OPPOSITE  IN  DIRECTION,  IN- 

C  CREASED  IN  LENGTH  EY  THE  FACTOR  1.3,  AND  COLLINEAB  WITH 

C  THE  REPLACED  VERTEX  AND  CENTEOID  OF  ALL  OTHER  VERTICES.) 

C 

C  WHEN  AN  OVER-REFLECTION  IS  NOT  FEASIBLE  OR  REMAINS 

C  WORST,  IT  IS  CONSIEERSD  NOT-PERMISSIBLE 

C  AND  IS  DISPLACED  EALFWAY  TOWARD  THE  CENTROID. 

C  AFTER  FOUR  SUCH  ATTEMPTS  ARE  MADE  UNSUCCESSFULLY 

C  EVERY  FIFTH  ATTEMPT  IS  MADE  BY  REFLECTING  THE  OFFENEING 

c  VERTEX  THROUGH  THE  PRESENT  EEST 

C  VERTEX,  INSTEAD  OF  THROUGH  THE  CENTROID.  IF  5*n+10 

C  DISPLACEMENTS  AND  CVER-REFLECTIONS  OCCUR  WITHOUT  A 

C  SUCCESSEUL  (PERMISSIBLE)  RESULT,  THE  CURRENT  BEST 

C  VERTEX  IS  TAKEN  AS  AN  INITIAL  FEASIBLE  POINT  FOR  A 

c  RESTABT  RUN  OF  THE  COMPLETE  PROCESS. 

C  RESTARTING  IS  ALSO  UNDERTAKEN  WHEN  6*nv+10  CONSECUTIVE 

C  TRIALS  HAVE  BEEN  MADE  WITH  NO  SIGNIFICANT  CHANGE  IN  THE 

C  VALUE  OF  THE  OBJECTIVE  FUNCTION.   IN  ALL  CASES, 

C  RESTARTING  IS  INHIEITED  IF  THE  LAST  RESTART  DID  NOT 

C  PRODUCE  A  SIGNIFICANT  IMPROVEMENT  IN  THE  MINIMUM 

C  ATTAINED. 

C 
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C    IT    IS    EECOMHENDED    THAT    THE    USEE    BEAD    THE    BEFEBENCE    FOB 
C    FDBTHEB    DSEFUL    INFOfiflATIOH-       IT    SHOULD    BE    NOTED    THAT    THE 
C    ALGORITHM    DEFINED    THEBE    HAS    BEEN    ALTERED    TO    FIND    THE 
C   CCNSTEAINED    MINIMUM,    BATHEB    THAN    THE    MAXIMUM. 
C 

c 
c 

C  REMARKS 

C 

C  THE  INTEGER  PROGRAMMING  OPTION  HAS  ADDED  TO  THIS  PROGRAM 

C  AS  SUGGESTED  IN  REFERENCE  (2)  .   A  MIXED 

c  integer/continuous  variable  version  of  boxplx 

C  WOULD  BE  EASY  TO  CREATE  BY  DEclaring  »ip»  to  be  an  array 

C  OF  NV  CONTROL  VARIABLES  WHERE  IP  (i) = 1  would  indicate 

C  THAT  THE  I-TH  VARIABLE  IS  TC  BE  CONFINED  to  integer 

C  VALUES.   EACH  STATEMENT  OF  THE  FOBM  'IF  (IP  . EQ- T) ■  etc. 

C  WOULD  THEN  NEED  TC  EE  ALTERED  TO  ■ IF  (IP  (I)  - EQ.  1)  •  etc. 

C  ,  WHERE  THE  SUBSCRIPT  IS  APPROPRI ATELY  CHOSEN.   NORMALLY, 

C  XU  AND  XL  VALUES  ARE  ALTERED  TO  BE  AN  EPSILON  'WITHIN* 

c  actual  values 

C  DECLARED  BY  THE  USER.   THIS  ADJUSTMENT  IS  NOT  MADE 

C  WHEN  IP=1. 

C 

C  NOTE:   NO  NON-LINEAR  PBOGRAMMING  ALGORITHM  CAN  GUARANTEE 

c  that  the  answer  found  is  the  global 

C  MINIMUM,  RATHER  THAN  JUST  A  local  minimum,  however, 

C  ACCORDING  TO  REF.2,  THE  COMPLEX  method  has  an  advantage 

C  IN  THAT  IT  TENDS  TO  FIND  THE  GLOBAL  minimum  more 

C  FREQUENTLY  THAN  MANY  OTHER  NCN-LINEAB  PROGRAM- 

C  MING  ALGORITHMS. 

C 

C  IT  SHCULD  BE  NOTEI  THAT  THE  AUXILIARY  VABIABLE  FEATURE 

c  CAN  AISO  BE  USED  1C  DEAL  WITH 

C  PRCELEMS  CONTAINING  EQUALITY  CONstraints.  any  equality 

C  CONSTRAINT  IMPLIES  THAT  A  GIVEN  VARiable  is  not  truly 

C  INDEPENDENT.   THEBEFOBE,  IN  GENEBAL,  ONE  variable 

C  INVOLVED  IN  AN  EQUALITY  CONSTRAINT  CAN  BE  RENUMBERED  from 

C  THE  SET  OF  NV  INDEPENDENT  VABIA3LES  AND  ADDED  TO  THE  SET 

C  OF  NAV  AUXILIABY  VABIABLES.   THIS  USUALLY  INVOLVES 

C  renumbering  THE  INIEPENDENT  VARIABLES  OF  THE  GIVEN 

C   problem 

C    SUBROUTINES    AND    FUNCTIONS    BEQUIBED 

C 

C  SUBROUTINE  'BOUT'  AND  FUNCTION  'FBV1  ABE  INTEGRAL 

C  parts  of  THE  BOXPLX  PACKAGE. 

C 

C  THO  FUNCTIONS  MUST  BE  SUPPLIED  BY  THE  USEE.   THE  FIRST, 

c  ke  (x)  .  is  used  to  evaluate  the  implicit 

C  CONSTRAINTS.   SET  KE=0  AT  THE  beginning  of  the  function 

C  ,  TEEN  EVALUATE  TEE  IMPLICIT  CONstraints.  in  the  example 

C  AECVE,  THE  FIEST  CONSTRAINT,  X (3)  ,  must  be  within  the 

C  EANGE  {0.  .LE.  X(3)  .LE.  6.).   THE  SECOND  constraint  x  (4) 

C  ,  MUST  BE  .GE.  0.  .   IF  EITHER  CONSTRAINT  IS  not  within 

C  THESE  BOUNDS,  CONTROL  IS  TBANSFEEBED  TO  STATEMENT  1. 

C  ANE  KE  IS  SET  TO  "1"  AND  CONTBOL  IS    BETURNED  TO  BOXPIX. 

C 

C  THE  SECOND  FUNCTION  THE  USER  MUST  PROVIDE  EVAIUATES  THE 

c  objective  function,  it  is 

C  CALLED  FEJX)  AS  SHOWN  IN  THE  EXAMple  above,  and  fe 

C  MUST  EE  SET  TC  THE  VALUE  OF  THE  OBJECTIVE  function 

C  CCBBESPCNDING  TO  CUEBENT  VALUES  OF  THE  NV  INDEPENDENT 

C  VABIAELES  IN  ARBAY  'X'. 

C 

C      BEFEBENCES 

C 

C  BOX,  M.  J.,  "A  NEW  METHOD  OF  CONSTBAINED  OPTIMIZATION 

C  and  a  COMPARISON  WITH  OTHER  METHODS", 

C  computer  journal,  8  apr.  »65,  PP.  45-52. 
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C  BEVERIDGE  G.,  AND  SCHECHTSR  R.  ,  "OPTIMIZATION:  THEORY  AND 

C  PRACTICE",  MCGRAW-EILL,  1970. 

C 

C  EBGGEAMMEfi 

C 

C  R.R.  HILLEAEY  1/1966. 

C  REVISED  TOR  SYSTEM  360  4/1967 

f*  CflRRFCTFn   1/1969 

C         REVISED/EXTENDED  BY  L. NOLAN/R. HILLEARY   2/1975 

C         CORRECTED   8/1976 

C 

c 

c 
c 

SUBROUTINE  BOXPLX  (NV,  NAV , NPR, NT2 ,RZ, XS , IP, BU ,BL , YMN ,IER) 
C 

DIMENSION  7(50,50),  FUN  (50) ,  SUM  (25),  CSN  (25) ,  XS(NV) 
1,tu  (nv) ,bl  (NV) 
C 

KV  =  5 

EE  =  1.E-6 

NTA  =  2000 

IF  (NTZ.GT.0)  MA  =  NTZ 

fl  =  RZ 

IF  (R.LE.0..OR.R.GE.  1.)  R=1./3. 

NVT  =  NV+NAV 
C 
C  TOTAL  VSfiS,  EXPLICIT  PLUS  IMPLICIT 

NT  =  0 
C  CURRENT  TRIAL  NO- 

NET =  0 
C  CURRENT  NO.  OF  PERMISSIBLE  TRIALS 

NTFS  =  0 
C         CURRENT  NO.  OF  TIMES  F  HAS  BEEN  ALMOST  UNCHANGED 
C 

C  CHECK  EEASIBILITY  OF  START  POINT 

C 

DO  4  1=1, NV 

VI  =  XSJI) 

IE  _(BL(I)  .LE.VT)  GO  TO  1 

VI  =  BL  (I) 
GC  TO  2 

1  IE  (BU(I)  .GE-VT)  GO  TO  3 
II  =  I 

VI  =  BU  (I) 

2  IF  (NPR.GT.0)  HRITE  (6,49)  II 

3  YJI,1)  = 


c 
c 


VJI,1)  =  VT 

CEN  (I)  =  VT 

IE  JIP.EQ.  1)  G< 

BI(I)  =  Bill)  +AMAX1  (EP,EP*ABS  (BL  II)  )  ) 

BD(l(  =  BUJI) -AMAX1 (EP,EP*ABS  (BU  (I) ) ) 

SUM(I)  =  VT 


NCE  =  1 
C  NUMBER  CF  CONSTRAINT  EVALUATIONS 

I  =  1 

IF  (KE(V  (1,  1)  ).EQ.0)  GO  TO  5 

IF  JnPR.LE-O)  GO  TO  12 

WRITE  (6,50) 

GC  TO  12 
5  NEE  =  1 
C 
C    KUKEEE  OF  VERTICES  (K)  =  2  TIMES  NO.  OF  VARIABLES. 

K  =  2*NV 
C 
C    NUMBER  OF  DISPLACEMENTS  ALLOWED. 
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SUM  =  5*HV+10 
C 

C    NUMEEB  OF  CONSECUTIVE  TRIALS  WITH  UNCHANGED  FE  TO 
c    terninate. 

NCT  =  NLIM+NV 

AIPHA  =  1.3 

FF  =  K 

FKM  =  FK-1. 

BETA  =  ALPHA+  1. 
C 
C    INSDBE  SEED  OF  EANDOM  NUMEEB  GENERATOH  IS  ODD. 

IQB  =  B*1.E7 

IF  (MOD  (IQR,2).EQ.0)  ICE=IQR+101 
C 
C  SET  DP  INITIAL  VERTICES 

FDN(1)  =  FE(V(1,1)) 

YEN  =  FUN(1) 

6  FI  =  1. 
FUNCLD  =  FUN  (1) 

C 

DC    15  1=2, K 
FI  =  FI+1. 
LIMT  =  0 

7  LIMT  =  LIMT+1 
C 

C  END    CALCULATION    IF    FEASIEIE    CENTROID    CANNOT    BE    FOUND. 

IE     (IIMT.GE.  Nllfl)1     GO    TO    11 
C 

DC    8    J=1,NV 
C 
C         BANECM    NUMBER    GOEEATOE       (RANDU) 

IQR    =    IQR*65539 

IF     (IQR.LT.O)     IQR    =    IQR+2 147483647+ 1 

EQX  =  IQR 

££X  =  RQX*.4656613E-9 

VJJ,I)  =  BL  (J)+RQX*(BU(J)-BL(J)) 

IF  (IP.  EQ.1)  V  (J,I)=AINT  (V  (J,  I)  +.5) 

8  CCNTINUE 
C 

DO  10  L=1,NLIM 

NCE  =  NCE+1 

II  (KE(V(1,I)  )  .EQ.O)  GO  TO  13 


C 

c 


DO  9  J=1,NV 

VT  =  .5*(VUI 

IF  (IP.EQ.1j  M    =  AINT(VT*.5) 

S  CCNTINUE 

10  CCNTINUE 


11  IF  (NPR.LE.O)  GO  TO  12 
WBITE  (6,51)  I 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,I,FUN,CEN,I) 

12  IER  =  -1 
GC  TO  48 

C 

13  DO  14  J=1,NV 

SUM  J)  =  SUM  (J)+V (J,I) 

14  CEN(J)  =  SUM  (0)/FI 

C    TRY  TO  ASSUEE  FEASIBLE  CENTROID  FOR  STARTING. 
NCE  =  NCE-H 

IF  (KE(CEN)  .EC..0)  GO  TO  60 
SUMJJ)  =  SUM  (J)  -V(J,I) 
GC  TO  7 
6C  NFE  =  NFE+1 
^c    FDNJIi  =  FE(V  (1,1)) 

15  CCNTINUE 
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c 

C  END    CF    LOOP    SETTING    OF    INITIAL    COMPLEX. 


IF    JNPR.LE.O)     GO    TO    17 

;nt,- 


CALL    BOUT     (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 

c 

C         FINE    THE    WORST    VIETEX,    TEF    »J'TH. 

J    =    1 
C 

DC    16   I  =  2,K 

IF_  (FDN  (J).GE.IUN  (I)  )     GC   TO    16 

16  CCNTINDE 
C 

C         EASIC    LOOP.       ELIMINATE    EACH    WORST    VERTEX    IN    TURN. 

C         it    must   become    NC    LONGER    WORST,    NOT    MERELY    IMPROVED. 

C         find    next-to-vertex,    THE    'JN'TH   ONE. 

17  JN    =    1 

IF     (J.EQ.  1)    JN   =    2 
C 

DC    18    1=1, K 

IF     (I-EQ-J)     GC    TO    18 

IF     (FUN(JN)  .GE.FUN(I)  )     GO    TO     18 

JN    =    I 

18  CONTINUE 
C 

C         1IMT    -    NUMBER    OF    MOVES    DURING    THIS    TRIAL   TOWARD    THE 
C  centxoid    DUE    TO   FUNCTION    VALUE. 

LIMT   =    1 
C 

C    COMPUTE  CENTROII  AND  OVEE  REFLECT  WORST  VERTEX. 
C 

DC  19  1=1, NV 

VI  =  V(I,  J) 

=  SUM 


SUM  (I)  =  SUM(I)-VT 

CEN  ij  =  SUM  (ij/FKM 

VT  =  BETA*CEN  (I)  -ALPHA; 


{I)  -ALPHA*VT 
IF  (IP.EQ.1)  \I  =  AINT(VT+.5) 
C 
C    INSURE  THE  EXPLICIT  CONSTRAINTS  ARE  OBSERVED. 

19  V(I,J)  =  AHAX1  (AMIN1  (VT,BU  (I)  )  ,BL  (I)  ) 
C 

NT  =  NT+1 
C 

C    CHECK  FOR  IMPLICIT  CONSTRAINT  VIOLATION. 
C 

20  DC  25  N=1,NLIM 
NCE  =  NCE  +  1 

IF  (KE(V  {1,  J)  )  .EQ.O)  GO  TO  26 
C 

C    EVERY  »KVfTH  TIME.  OVER-REFLECT  THE  OFFENDING  VERTEX 
C    through  the  BESI  VERTEX. 

IF  (MOD  (N,KV)  -KE.0)  GO  TO  22 

CALL  FBV  (K,F0fi,H) 
C 

DC  21  1=1, NV 

VT  =  BETA*V  (I.E)-ALPHA*V(I,J) 

IF  (IP.EQ.1)  \I  =  AINTJVT+.5) 

21  V(I,J)  =  AMAX1  (AMIN1  (VT,BU{I)  )  ,BL  (I)) 
C 

GC  TO  24 
C 

C    CCNSTRAINT  VIOLATION:   MOVE  NEW  POINT  TOWARD  CENTROID. 
C 


22  DC  23  1=1 ,NV 

VT  =  .5*  (CEN  (I) +V  (I,  J)) 

IF  (IP.EQ.  1)  Ml   =  AINT( 


VJI,J)  =  VT 

""NU" 


(IP.EQ.1)  M  =  AINT(VT+.5) 
23  CONTINUE 
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c 
C 
C 


C 

c 


c 
c 


c 
c 
c 
c 


c 
c 

c 
c 
c 
c 


c 
c 


c 
c 
c 

c 


24  NT  =  NT+1 

25  CONTINUE 

I1B  =  1 

CANNOT  GET  FEASIELE  VERTEX  BY  MOVING  TOWARD  CENTROID, 
CB  EY  OVER-REFLICTING  THRU  THE  BEST  VERTEX. 

IF  (NPR.LE-0)  GO  TO  42 

WBITE  (6,52)  NT, J 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,J) 

GC  TO  42 

FEASIELE  VERTEX  FOUND,  EVALUATE  THE  OBJECTIVE  FUNCTION 

26  NFE  =  NFE  +  1 

FUNTRY  =  FE{V  (1,J))  "» 

TEST  TO  SEE  IF  FUNCTION  VALUE  HAS  NOT  CHANGED, 
AFO  =  ABS  (FUNTEY-FUNOLD) 
AMX  =  AMAX1  (AES (EP*FUNOLD)  ,EP) 

ACTIVATE  THE  FOLIOWING  TWO  STATEMENTS  FOR  DIAGNOSTIC 
purposes  only. 

WRITE  (6,99)  J, AFO, AMX, FUNTRY, FUNOLD, FUN (J) ,FUN  (JN) 
1,NTFS,N 
99  FORMAT  ( 1  X,  13  ,6E1  5.7,  215) 

IF  (AFO. GT. AMX)  GO  TO  27 

NTFS  =  NTFS+1 

IF  (NTFS.LT.  NCT)  GO  TO  28 

IFR  =  0 

IF  (NPR.LE.O)  GO  TO  42 

WRITE  (6,53)  K 

CALL  BOUT  (NT,NPT,NFE,NCE,NV, NVT, V, K, FUN,CEN, 0) 

GO  TO  4  2 

27  NTFS  =  0 

IS  TEE  NEW  VERTEX  NO  LCNGFE  WORST? 

28  IF  (FUNTRY. LT. FUN  (JN)  )  GO  TO  34 

TRIAL  VERTEX  IS  STILL  WORST:  ADJUST  TOWARD  CENTEOIB. 
EVERY  'KV»TH  THE,  OVER-REFLECT  THE  OFFENDING  VERTEX 
through  the  BEST  VERTEX. 

LIMT  =  LIMT+1 

IF  (MOD  (LIMT,  KV)  .NE.  0)  GO  TO  30 

CALL  FBV  (K,FU*,M) 

DC  29  1=1, NV 

VT  =  BETA*V(I,M)  -ALPHA*V  (I, J) 

IF  (IP.  EQ.  1)  VI  =  AINTJVT  +  .5' 

29  V(I,J)  =  AMAX1  (AMIN1  (VT,BU  (I)  )  ,BL  (I)  ) 

GC  TO  32 

30  DO  31  1=1. NV 
VT  =  .5*  (CEN(I)+V  (I.J)) 
IF  (IP.EQ.  1)   \1  =  AINTf 


(IP.EQ.1)   \1  =  AINT(VT+.5) 
VJ1-J)  =  VT 

31  CONTINUE 

32  IF  (LIMT.LT.NLIM)  GO  TO  33 


CANNOT  MAKE  THE  »J'TH  VERTEX  NO  LONGER  WORST  BY 
displacing  toward 
THE  CENTROID  OR  BY  OVER-REFLECTING  THRU  THE  BEST  VERTEX. 
IER  =  2 

IF  (NPR  .LE.  0)   GO  TO  42 
WBITE  (6,52)   KT,  J 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 
GC  TO  42 
33  NT  =  NT+1 
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GC  TO  20 
C 

C    SUCCESS:   WE  HAVI  A  REPLACEMENT  FOE  VERTEX  J. 
34  FOB  (J)  =  FUNTBY 

FUNOLD  =  FUNTEI 

NPT  =  NPT+1 
C 

C    EVERY  100»TH  PERMISSIBLE  TRIAL,  RECOMPUTE  CENTRCID 
C    summation  to  AVCID  CREEPING  ERROR. 

IF  (MOD  (NPT,1CC)  .NS.  0)  GO  TO  37 


C 

c 
c 

c 


c 

c 


DO  36  1=1, NV 
SUM(I)  =  0. 

DC  35  N=1 ,K 

35  SUM  (I)  =  SUM  (I)  +  V(I,N) 

CEN  (I)  =  SUM  (I)/FK 

36  CCNTINUE 

LC  =  0 
GO  TO  3  9 

37  DO  38  1=1, NV 

38  SUM  (I)  =  SUM  (I)+V (I, J) 

LC  =  J 

39  IF  (NPR.LE.O)  GO  TO  40 

IF  (MOD  (NPT,NPB) -NE.0)  GO  TO  40 


CALL  BOUT  (NT,N?T,NFE,NCS,NV,NVT,V,K,FUN,CEN,LC) 
C 

C    EAS  THE  MAX-  NUMEER  OF  TRIALS  BEEN  REACHED  WITHOUT 
C    convergence?  iF  KOT,  GO  TO  NEW  TRIAL. 

40  IF  (NT.GE.NTA)  GO  TO  41 
C 

C    NEXT-TO-WORST  VERTEX  NOW  BECOMES  WORST. 
J  =  JN 
GC  TO  17 

41  IEB  =  3 

IF  (NPR.GT.0)  RRITE  (6,54) 
C 

C    COLIECTOR  POINT  FOR  ALL  ENDINGS. 

C   1)   CANNOT  DEVELCE  FEASIBLE  VERTEX.  IER  =  1 

C  2)  CANNOT  DEVELOE  A  NO- LCNGER-WORST  VERTEX.  IER  =  2 
C  3}  FUNCTION  VALUE  UNCHANGED  FOR  K  TRIALS.  IER  =  0 
C   4]   LIMIT  ON  TRIALS  REACHED.  IER  =  3 

C   51   CANNOT  FIND  JEASISLE  VERTEX  AT  START.      IER  =  -1 

42  CCNTINUE 
C 

C    FIND  BEST  VERTEX. 

CALL  FBV  (K.FU1<,M) 

IF  (IER.GE.3)  GO  TO  44 
C 

C    RESTART  IF  THIS  SOLUTION  IS  SIGNIFICANTLY  BETTER  THAN 
C    the  previous,  OR  IF  THIS  IS  THE  FIRST  TRY. 

IF  JNPR.LE.0)  GO  TO  43 


WRITE  (6,55)  (P.YMN.FUN  (M)  ) 
IF  (FUN  (M).GE-YflN)  GO  TO  47 
IF  (ABS  (FUN (M)-YMN).LE.AMAX1 


43 

(EP,EP*YMN))  GO  TO  47 
C 
C    GIVE  IT  ANOTHER  TRY  UNLESS  LIMIT  ON  TRIALS  REACHED. 

44  YMN  =  FUN(M) 


FUN(1)  =  FUN  (M) 

DC  45  1=1, NV 
CEN  (I)  =  V(I,M) 
SUM(I)  =  V(I,M) 
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c 
c 


45  V(I,1)  =  V(I,K) 

DC  46  1=1  ,NVT 

46  XS(I)  =  V  (I,M) 

If  (IER.LT.3)  GO  TO  6 

47  IF  (NPR.LE.O).  GO  TO  48 
CALL  BOUT  (NT,NPT,NFE, 


(NT,NPT,NFE/NCE,NV,NVT,V, K,FUN,V  (1,M)  ,-1) 


WRITE  (6,56)  FUN(M) 
48  RETURN 

49  FORMAT  (50H0INDEX  AND  DIRECTION  OF  OUTLYING 
1  variable  at  starti5) 

50  FORMAT  (50  HO  IMPLICIT  CONSTRAINT  VIOLATED  AT 
Istart.  dead  end.) 

51  FCEMAT  ('OCANNCT  FIND  FEASIBLE' , 14, »TH  VERTEX  OR 
Icentroia  at  start.'} 

52  FORMAT  (10H0AI  TRIAL  I4,54H  CANNOT  FIND  FSASIELE 
Ivertex  which  is  no 

2LCNGER  WORST, 14, 15X, 'RESTART  FROM  BEST  VERTEX.') 

53  FORMAT  J4 0H0F0KCTION  HAS  BEEN  ALMOST  UNCHANGEE 
1for  i5,7h  trails) 

54  FCRKAT  {27H0LIMT  ON  TRIALS  EXCEEDED.  ) 

55  FCRMAT('03EST  VERTEX  IS  NO. ',13, 'OLD  MIN  WAS',E15.7, 
1  '   NEW  MIN  X5  '  E15.7) 

56  FORMAT  (• OMIN  OBJECTIVE  FUNCTION  IS  ',E15.7) 
END 

SUEROUTINE  FBV  (K,FUN,M) 
DIMENSION  FUN  (50) 

a  =  1 

DC  1  1=2, K 

IF_  (FUN(M).LE.FUN(I)  )  GOTO  1 

1  CONTINUE 

RETURN 

END 

SUBROUTINE  BOUT  (NT, NPT , NFE, NCE, NV, NVT,V,K ,FN ,C, IK) 

DIMENSION  V(50,50),  FN  (50),  C (25) 

WRITE  (6,4)  NI,NPT,NFE,NCE 

DC  1  1=  1  K 

WRITE  (6*5)  FN  (I)  ,  (V{J,I)  ,J=1,NV) 

IF  (NVT.LE.NV)  GO  TO  1 

NVP  =  NV+1 

WHITE  (6,6)   (V(J,I)  ,J=NVE,NVT) 

1  CONTINUE 

IF  (IK.NE.O)  GO  TO  2 

WRITS  (6,7)   (C(I)  ,I=1,NV) 
RETURN 

2  IF  (IK.GE.O)  GC  TO  3 
WRITE  (6,8)   (C(I)  ,I=1,NV) 
EFTURN 

3  WRITE  (6,9)  IK,  (C  (I)  ,1=1, NV) 
BETURN 

4  FORMAT  (»0NO.  TOTAL  TRIALS  =  »,I5,4X, 
1'no.  feasible  trails  =  ',i5,4x, 
2'NO.  FUNCTION  EVALUATIONS  =  »,I5,4X, 
3'nc.  constraint  evaluations  =  ',i5/ 

4»0      FUNCTION  VALUE ' ,6X  'INDEPENDENT  VARIABLES/ 
5dependENT  OR  IMPLICIT  CONSTRAINTS') 

5  FORMAT     (1H    ,E18.7,2X,7E14.7/(21X,7E14.7)) 

6  FORMAT     (21X.7E14.7) 

7  FORMAT     (10H0CFNTROID    1 1 X, 7E1 4 .7/ (21 X. 7E14.7) ) 

8  FORMAT     («0       BEST    VERTEX » ,7X, 7E14.7/ (21X, 7E14. 7) ) 
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9  FCEMAT  ('OCENTEOID  LESS  VX' , 12 ,2X, 7E 14.7/ (2 1 X, 
17e14.7)  ) 
END 

FUNCTION  FE{X) 
DIMENSION  X{3) 
CCHMCN  TDIFF 
CALL  PLANT  (X) 
FE=TDIFF 
BETURH 
END 

FUNCTION  KE(X) 
DIMENSION  X{3) 
KE=0 
BETUEN 
END 
//GO.SYSIN  DD  * 

/* 

//GO.ET12F001  DD  DISf=SHH#DSN-MSS. S2160. A341 


85 


APPENDIX  B 
EXAflPIE  PROBLEM  OSIHG  ICSOS 

The  purpose  of  this  example  is  to  demonstrate  the 
performance  of  the  program.  Consider  the  control  syst€m  of 
Figure  4.1  with  controller  C.  Figure  B.  1  shows  the  fclock 
diagram  to  evaluate  the  controller  parameters. 
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Ihe  differential  equations  describing  the  system  and  its 
desired  performance  are: 

22=DX2 

X4=DX4 

XS=DX5 

fi  =DR 

YAH=R 
Defining  the  following  cost  function: 

J=  J (LAMDA*YAHi**2+D**2) dt 
Defining  the  special  functions: 

YAHE=YASC-YAW 

EX2=(YAWE-X2)  /12 

X3=K1*{T1*DX2+X2) 

DX4=  (X3-X4J/T4 

d=(t3*dx4+x4) 

dx5=  (d-x5)/tp1 

x8=k*(tz*dx5  +  x5) 

dr=  (x8-r)  /tp2 
Defining  the  constants: 

YAWC=1.0 

K=.  14771 

12=11.2833 

1P1=6.4699 

1£2=53.7931 

LAMDA=4.2 
and  using  YAWC=1.0  the  optimal  solution  found  by 
the  program  is: 

K1=. 4179916 

11=53.69932 

12=4.970023 

13=6.294369 

14=13.85919 

CCS!=68. 04735 
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Table  34  shows  the  specifications  of  this  problem  with 
the  free  parameter  optimum  values  found.  Figure  B.2  shows 
the  actual  yaw  and  rudder  response. 
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TABLE   34 
ICSOS    ODTPDT 

SPECIFIC AT ICNS 

VARIABLES    £     INITIAL    CCNCITIONS: 

X2    =    .0 

X4    =    ,0 

X5    =    .0 

R    =    .0 

YAW    =    .C 

J    =    .0 

TIME   =    .0 

CONSTANTS: 
YAWC    =    l.OOOOOCOOO 
K    =    .147710000C 
TZ    =    11.233300CC 
TP1    =    6.469930COO 
TP2    =    52o793100C0 
LAMDA    =    4.2OC0CC0O0 


FREE    PARAMETERS: 

Kl    : 

:    CV=    .4179916 

LL  = 

.lcccooo 

UL  = 

1.00000  0 

Tl    : 

:    0V=    53. 69*28 

LL  = 

1*000000 

LL  = 

100,0000 

T2    ! 

:    0V=    4.970C29 

LL  = 

l.QOCOOO 

UL= 

!C.  0000  0 

T3    : 

:    0V=    6.294369 

LL  = 

1.000000 

UL  = 

10.0000  0 

T4    : 

I    0V=    13.85917 

LL  = 

l.COCOOO 

LL  = 

2C.0OOOC 

SPECIAL    FUNCTICNS: 
YAWE    =    YAWC-YAV» 
DX2    =    (YAWE-X2)/T2 
X3    =    K1*(T1*DX24X2 ) 
0X4    =    (X2-X4)/T4 
D   =     (T3*CX4+X4 ) 
DX5    =    (C-X5) /TF1 
X8    =    K*(TZ*DX5+X5) 
OR    =    U8-R)/TP2 

DERIVATIVES: 

D<X2    /D(TIME     )    =   = 

0X2 
0(X4    /DCTIME     )    =    = 

DX4  N 

DJX5    /DUIME     )    =    = 

DX5 
D<R    /DiTIME     )    =   = 

OR 
D(YAK    /C(TIME     )    =    = 

R 
D(J    /0(TIME    )    =   = 

LAM0A*YAWE**2+D**2 

OUTPUTS: 

TITLE:    ACTUAL   YAW    ANC    RUDDER    RESFCNSE 
TABULATE:    TINE  D  R  YAW 

AT    INTERVAL  2.CC0COO000 

PLOT:    D  YAW 

AGAINST:    TIME         AT    INTERVAL  2.000000000 

END    CALCULATION   WHEN    TINE    oGE.     600. OCO 
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